Index: /trunk/tools/DistriCorrection.f90
===================================================================
--- /trunk/tools/DistriCorrection.f90	(revision 1609)
+++ /trunk/tools/DistriCorrection.f90	(revision 1609)
@@ -0,0 +1,94 @@
+PROGRAM DistriCorrection
+! Empty program to compuile the modules for the modification of a Distribution
+! L. Fita, CIMA August 2017
+!
+!! This is because I do not know enough about 'Makefile'
+
+  USE module_generic
+  USE module_scientific
+  USE module_DistriCorrection
+
+  IMPLICIT NONE
+
+  INTEGER                                                :: dim1, dim2, dim3, dimrng, dimpercen, Numquants
+  REAL(r_k)                                              :: Highquant, fValue
+  REAL(r_k), DIMENSION(:,:,:), ALLOCATABLE               :: values
+  REAL(r_k), DIMENSION(:), ALLOCATABLE                   :: rangeValues
+  REAL(r_k), DIMENSION(:), ALLOCATABLE                   :: PercentagesRange
+  CHARACTER(len=50)                                      :: KRemove, ReplacevalS
+  LOGICAL                                                :: debug
+  REAL(r_k), DIMENSION(:,:,:), ALLOCATABLE               :: newvalues
+  REAL(r_k), DIMENSION(:,:,:), ALLOCATABLE               :: rangeamount
+  INTEGER, DIMENSION(:,:,:), ALLOCATABLE                 :: redistributedindices, aboveindices
+  REAL(r_k), DIMENSION(:,:), ALLOCATABLE                 :: Highabovequant
+
+! Local
+  INTEGER                                                :: i,j
+
+!!!!!!! Variables
+! dim[1/2/3]: size of the values to treat
+! values: values to re-distribute
+! fValue: value for the missing ones
+! dimrng: number of ranges provided
+! rangeValues: ',' values from which establish the ranges of removals
+! dimpercen: number of percentages provided
+! PercentagesRange: ',' percentages of removing to apply to each range (one value less than [rangeVals])
+! KRemove: kind of methodology to select which values to remove
+!   'rand': randomly selection within the correspondent range
+! Numquants: quantity of percentiles to compute
+! Highquant: percentile of the distribution of values above which total removed amount will be equi-distributed
+! ReplacevalS: replacing value to be assigned at the removed values
+!   'min': initial minimal value of the distribution
+!   'mean': initial mean value of the distribution
+!   '[value]': any given value otherwise
+!
+! Outcome:
+! newvalues: Values of the new distribution
+! redistributedindices: List of positions changed for each range
+! rangeamount: Amount re-distributed from each range
+! Highabovequant: value from which the amount is distributed
+! aboveindices: Indices of the values which received the distributed amount
+
+  fname = 'DistriCorrection'
+
+  dim1 = 5
+  dim2 = 3
+  dim3 = 100
+  dimrng = 3
+  dimpercen = 2
+  fValue = -99999.*oneRK
+
+  ALLOCATE(values(dim1,dim2,dim3))
+  ALLOCATE(rangeValues(dimrng))
+  ALLOCATE(percentagesRange(dimpercen))
+  ALLOCATE(newvalues(dim1,dim2,dim3))
+  ALLOCATE(rangeamount(dim1,dim2,dimrng))
+  ALLOCATE(redistributedindices(dim1,dim2,dim3))
+  ALLOCATE(aboveindices(dim1,dim2,dim3))
+  ALLOCATE(Highabovequant(dim1,dim2))
+  
+  DO i=1,dim1
+    DO j=1,dim2
+      values(i,j,:) = RangeR_K(dim3, oneRK, dim3*oneRK)
+    END DO
+  END DO
+  PRINT *,values(1,1,:)
+
+  rangeValues = (/ 0., 5., 25./)
+  rangeValues = (/ 0., 5., 25./)
+  percentagesRange = (/ 0.8, 0.5/)
+
+  KRemove = 'rand'
+  Numquants = 20
+  Highquant = 0.9d0
+  ReplaceValS = '0.'
+  debug = .TRUE.
+
+  ! Testing the module
+
+  CALL randLowRemovePercent_HighEqui3D(dim1, dim2, dim3, values, dimrng, rangeValues, dimpercen,      &
+    PercentagesRange, KRemove, Numquants, Highquant, ReplacevalS, fvalue, debug, newvalues,           &
+    redistributedindices, rangeamount, Highabovequant, aboveindices)
+
+END PROGRAM DistriCorrection
+
Index: /trunk/tools/module_DistriCorrection.f90
===================================================================
--- /trunk/tools/module_DistriCorrection.f90	(revision 1609)
+++ /trunk/tools/module_DistriCorrection.f90	(revision 1609)
@@ -0,0 +1,428 @@
+MODULE module_DistriCorrection
+! Modules for modIFication of a Ditribution
+
+  USE module_definitions
+  USE module_generic
+  USE module_scientific
+
+CONTAINS
+
+!!!!!!! Subroutines
+! randLowRemovePercent_HighEqui: Method to re-distribute from lower tail of thedistribution to the 
+!   higher one following percentages preserving total amount. 
+
+SUBROUTINE randLowRemovePercent_HighEqui3D(d1, d2, d3, dist, drng, rangeVals, dpercen, PercentRange,  &
+  KindRemove, Nquants, Hquant, RepvalS, fillValue, dbg, newdist, redistributed, rngamount,            &
+  Habovequant, aboveind)
+! Method to re-distribute from lower tail of thedistribution to the higher one following percentages 
+!   preserving total amount. This method, removes values from the lower end of the distribution and 
+!   equi-distribute the total amount on a percentile on the high end
+
+  IMPLICIT NONE
+
+  INTEGER, INTENT(in)                                    :: d1, d2, d3, drng, dpercen, Nquants
+  REAL(r_k), INTENT(in)                                  :: Hquant, fillValue
+  REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)             :: dist
+  REAL(r_k), DIMENSION(drng), INTENT(in)                 :: rangeVals
+  REAL(r_k), DIMENSION(dpercen), INTENT(in)              :: PercentRange
+  CHARACTER(len=50), INTENT(in)                          :: KindRemove, RepvalS
+  LOGICAL, INTENT(in)                                    :: dbg
+  REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out)            :: newdist
+  REAL(r_k), DIMENSION(d1,d2,drng), INTENT(out)          :: rngamount
+  INTEGER, DIMENSION(d1,d2,d3), INTENT(out)              :: redistributed, aboveind
+  REAL(r_k), DIMENSION(d1,d2), INTENT(out)               :: Habovequant
+
+! Local
+  INTEGER                                                :: i, j, iv, ir, irm, ind, ia
+  INTEGER                                                :: ierr
+  INTEGER                                                :: irg, erg
+  INTEGER                                                :: iHquant
+  INTEGER                                                :: Nvals, Nrngvals, Nrmvals, Nredis, Nrmprev
+  INTEGER                                                :: NaboveHquants
+  INTEGER, DIMENSION(d3)                                 :: distind
+  REAL(r_k), DIMENSION(d3)                               :: varvals
+  REAL(r_k), DIMENSION(:),ALLOCATABLE                    :: varvalues, quants, newquants
+  REAL(r_k)                                              :: one, zero
+  REAL(r_k)                                              :: ndist, xdist, meandist
+  REAL(r_k)                                              :: Repval, leakvals
+  REAL(r_k)                                              :: iamount, initial_amount, eamount
+  REAL(r_k)                                              :: TOTamount, dTOTamount
+  INTEGER, DIMENSION(:), ALLOCATABLE                     :: rngind, rmind
+  CHARACTER(len=4)                                       :: IS1, IS2
+  CHARACTER(len=50)                                      :: spc, RS1
+  CHARACTER(len=50)                                      :: fname
+  CHARACTER(len=1000)                                    :: vectorS1, vectorS2, bspc
+  LOGICAL, DIMENSION(d3)                                 :: search1, search2, search
+  CHARACTER(len=256)                                     :: msg
+
+!!!!!!! Variables
+! d[1/2/3]: size of the values to treat
+! dist: values to re-distribute
+! drng: number of ranges provided
+! rangeVals: ',' values from which establish the ranges of removals
+! dpercen: number of percentages provided
+! PercentRange: ',' percentages of removing to apply to each range (one value less than [rangeVals])
+! KindRemove: kind of methodology to select which values to remove
+!   'rand': randomly selection within the correspondent range
+! Nquants: quantity of percentiles to compute
+! Hquant: percentile of the distribution of values above which total removed amount will be equi-distributed
+! RepvalS: replacing value to be assigned at the removed values
+!   'min': initial minimal value of the distribution
+!   'mean': initial mean value of the distribution
+!   '[value]': any given value otherwise
+! fillValue: missing value
+!
+! Outcome:
+! newdist: Values of the new distribution
+! redistributed: List of positions changed for each range
+! rngamount: Amount re-distributed from each range
+! Habovequant: value from which the amount is distributed
+! aboveind: Indices of the values which received the distributed amount
+
+  fname = 'randLowRemovePercent_HighEqui'
+
+  spc = ' '
+  bspc = Nstrings(spc,4)
+  one = 1.
+  zero = 0.
+
+  IF (dbg) PRINT *,"  Re-aranging distribution of ", d3," values following '"//TRIM(fname)//"' method"
+
+  ! Checking for the right amount of ranges and percentages
+  IF (dpercen /= drng - 1) THEN
+    WRITE(IS1,'(I4)')drng
+    WRITE(IS2,'(I4)')dpercen
+    vectorS1 = vectorR_KS(dpercen, PercentRange)
+    vectorS2 = vectorR_KS(drng, rangeVals)
+
+    msg = 'provided number of percentages: ' // TRIM(vectorS1) // ' (' // TRIM(IS2) //                &
+         ') does not fit for: ' // TRIM(IS1) // 'ranges:' // TRIM(vectorS2)
+    CALL ErrMsg(msg, fname, -1)
+  END IF
+
+  IF (dbg) THEN
+    PRINT *,TRIM(bspc) // 'Ranges modIFication _______'
+    DO ir=1, drng-1
+      WRITE(RS1,'(F5.2)')PercentRange(ir)*100
+      PRINT *,TRIM(bspc) // '  ' // TRIM(RS1) // ' % of (', rangeVals(ir), ',', rangeVals(ir+1), ']'
+    END DO
+  END IF
+
+  IF (ALLOCATED(quants)) DEALLOCATE(quants)
+  ALLOCATE(quants(Nquants+1), STAT=ierr)
+  msg = "Problem allocating 'quants'"
+  CALL ErrMsg(msg, fname, ierr)
+
+  ! Checking for the upper-quantile
+  quants = RangeR_K(Nquants+1, zeroRK, 100.*oneRK)
+
+  iHquant = Index1DArrayR_K(quants, Nquants+1, hquant*100.*oneRK)
+
+  IF (iHquant == -1) THEN
+    WRITE(RS1,'(F50.10)')hquant*100.
+    WRITE(IS1,'(I4)')Nquants
+    vectorS1 = vectorR_KS(Nquants+1, quants)
+    msg = 'provided high percentile ' // TRIM(RS1) // ' does not fit with ' // TRIM(IS1) //           &
+      ' quantiles: ' // TRIM(vectorS1)
+    CALL ErrMsg(msg, fname, iHquant)
+  END IF
+
+  IF (dbg) PRINT *,TRIM(bspc) // 'To be added above:', quants(iHquant), ' [', iHquant, '] percentile' 
+
+  ! getting the right replacing value
+  SELECT CASE (TRIM(repvalS))
+    CASE ('min')
+      repval = ndist
+    CASE ('mean')
+      repval = meandist
+    CASE DEFAULT
+      READ(RepvalS,'(F50.25)', IOSTAT=ierr)Repval
+      msg = "For not 'min' or 'mean', a value is required as RepValS" // CHAR(10) // " provided: '" //&
+        TRIM(RepvalS) // "'"
+      CALL ErrMsg(msg, fname, ierr)
+  END SELECT
+
+  ! Future new distribution
+  newdist = dist
+
+  ! Total amounts for each range
+  rngamount = zero
+
+  ! Distribution indices
+  distind = rangeI(d3,1,d3)
+
+  ! Initialization of indices
+  redistributed = -1
+  aboveind = -1
+
+  !! Looping on all dimensions
+  DO i=1, d1
+    DO j=1, d2
+
+      varvals = dist(i,j,:)
+      search = (varvals == fillValue)
+      Nvals = COUNT(.NOT.search)
+
+      doingnothing: IF (Nvals == 0) THEN
+        PRINT *,TRIM(infomsg)
+        PRINT *,'  ' // TRIM(fname) // ': all missed values for point:', i,', ',j, '!!'
+        newdist(i,j,:) = dist(i,j,:)
+        CYCLE
+
+      ELSE
+
+        IF (ALLOCATED(varvalues)) DEALLOCATE(varvalues)
+        ALLOCATE(varvalues(Nvals))
+
+        varvalues = PACK(varvals, MASK=.NOT.search)
+
+        ! Amount of previous indices removed
+        Nrmprev = 0
+
+        ndist = MINVAL(varvalues)
+        xdist = MAXVAL(varvalues)
+        meandist = SUM(varvalues)/(d3*one)
+
+        IF (xdist < rangeVals(drng)) THEN
+          PRINT *,TRIM(infomsg)
+          PRINT *,'  ' // TRIM(fname) // ': all values for point:', i,', ',j, ' are below ' //        &
+            'ranges upper-end:', rangeVals(drng), ' !!'
+          newdist(i,j,:) = dist(i,j,:)
+          CYCLE 
+        END IF
+
+        CALL quantilesR_K(d3, varvalues, Nquants+1, quants)
+
+        iamount = SUM(varvalues)
+        initial_amount = iamount + 0.
+    
+        IF (dbg) THEN
+          PRINT *,TRIM(bspc) // 'values _______'
+          PRINT *,TRIM(bspc) // '  initial total amount:', iamount
+          PRINT *,TRIM(bspc) // '  min:', ndist, 'max:', xdist, 'mean:', meandist
+          PRINT *,TRIM(bspc) // '  quantiles:'
+          CALL PrintQuantilesR_K(d3, varvalues, Nquants+1, quants, Nstrings(spc,6))
+        END IF
+
+        IF (dbg) THEN
+          PRINT *, TRIM(bspc) // 'Loooping ranges _______'
+        END IF
+        bspc = Nstrings(spc,6)
+
+        ! Getting total amount to remove and remove it from the value
+        TOTamount = zero
+
+        ranges: DO ir=1, drng-1
+          irg = ir
+          erg = ir + 1
+
+          IF (dbg) THEN
+            WRITE(RS1, '(F6.2)')PercentRange(ir)*100
+            PRINT *,TRIM(bspc) // '* re-aranging ' // TRIM(RS1) // ' % of interval (', rangeVals(irg), ', ',  &
+              rangeVals(erg), ']'
+          END IF
+  
+          ! Searching the indices
+          WHERE (varvalues > rangeVals(irg)) 
+            search1 = .TRUE.
+          ELSEWHERE
+            search1 = .FALSE.
+          END WHERE
+  
+          WHERE (varvalues <= rangeVals(erg)) 
+            search2 = .TRUE.
+          ELSEWHERE
+            search2 = .FALSE.
+          END WHERE
+
+          WHERE (search1 .AND. search2)
+            search = .TRUE.
+          ELSEWHERE
+            search = .FALSE.
+          END WHERE
+
+          ! Values within the range
+          Nrngvals = COUNT(search)          
+          IF (ALLOCATED(rngind)) DEALLOCATE(rngind)
+          ALLOCATE(rngind(Nrngvals), STAT=ierr)
+          msg = "Problem allocating 'rngind'"
+          CALL ErrMsg(msg, fname, ierr)
+
+          irm = 1
+          DO iv=1, d3
+            IF (search(iv)) THEN
+              rngind(irm) = iv
+              irm = irm + 1
+            END IF
+          END DO
+
+          ! Number of values to remove
+          Nrmvals = INT(Nrngvals*PercentRange(ir))
+
+          IF (Nrmvals < 2) THEN
+            PRINT *,TRIM(warnmsg)
+            PRINT *,'  ' // TRIM(fname) // ': very few values to remove:', Nrmvals, ' !!'
+            PRINT *,'    for '//TRiM(RS1)//' % of interval (', rangeVals(irg), ',', rangeVals(erg), ']'
+          END IF
+
+          ! Indices to remove
+          IF (ALLOCATED(rmind)) DEALLOCATE(rmind)
+          ALLOCATE(rmind(Nrmvals), STAT=ierr)
+          msg = "Problem allocating 'rmind'"
+          CALL ErrMsg(msg, fname, ierr)
+
+          ! Removing kind
+          Kremove: SELECT CASE (TRIM(KindRemove))
+  
+            CASE ('rand')
+              CALL rand_sample(Nrngvals, Nrmvals, rmind)
+
+            CASE DEFAULT
+              msg = "Removing kind '" // TRIM(KindRemove) // "' not ready" // CHAR(10) //               &
+                "Available ones: 'rand'"
+              CALL ErrMsg(msg, fname, -1)
+  
+          END SELECT Kremove
+
+          leakvals = zero
+          DO irm=1, Nrmvals
+            ind = rngind(rmind(irm))
+            ! keeping the indices of the redistributed data
+            redistributed(i,j,Nrmprev+irm) = ind
+            rngamount(i,j,ir) = rngamount(i,j,ir) + (varvalues(ind) - repval)
+            ! Keeping track of what is lost or gain during the transformation (for repval /= 0.)
+            leakvals = leakvals + (repval - varvalues(ind))
+
+            newdist(i,j,ind) = repval
+          END DO
+
+          TOTamount = TOTamount + rngamount(i,j,ir)
+          IF (dbg) THEN
+            PRINT *,TRIM(bspc) // '  removing:', Nrmvals, 'values from:', Nrngvals
+            PRINT *,TRIM(bspc) // '  from indices:', rngind, 're-distributing:', rngind(rmind(1:Nrmvals))
+            PRINT *,TRIM(bspc) // '  total amount to remove:', TOTamount
+          END IF
+          Nrmprev = Nrmprev + Nrmvals
+
+        END DO ranges
+
+        NonullRM: IF (TOTamount /= zero) THEN
+          ! Equi-distributing the TOTamount on values >= hquant
+          Habovequant(i,j) = quants(iHquant)
+  
+          IF (Habovequant(i,j) < rangeVals(drng)) THEN
+            PRINT *,TRIM(infomsg)
+            PRINT *,'    point:', i,', ',j,' has a High-quantile:', quants(iHquant), ' smaller' //    &
+              ' than upper limit of ranges:', rangeVals(drng), '!!'
+            !PRINT *,'    Trying to re-distribute on all the values above:', rangeVals(drng)
+            Habovequant(i,j) = oneRK
+
+            ! Searching the indices
+            !WHERE (varvalues >= rangeVals(drng)) 
+            WHERE (varvalues > oneRK) 
+              search = .TRUE.
+            ELSEWHERE
+              search = .FALSE.
+            END WHERE
+            PRINT *,'    Re-distribute ', TOTamount, ' on', COUNT(search),' values above:', oneRK
+
+          ELSE
+            ! Searching the indices
+            WHERE (varvalues >= Habovequant(i,j)) 
+              search = .TRUE.
+            ELSEWHERE
+              search = .FALSE.
+            END WHERE
+          END IF
+      
+          ! Quantity to add at each above quantile
+          NaboveHquants = COUNT(search)
+          dTOTamount = TOTamount / NaboveHquants
+  
+          IF (dbg) THEN 
+            PRINT *,TRIM(bspc) // 'incrementing:', NaboveHquants, 'values above:', Habovequant(i,j),  &
+              ' by', dTOTamount
+            PRINT *,TRIM(bspc) // '  replacing value:', repval
+          END IF
+
+          irm = 1
+          DO iv=1, d3
+            IF (search(iv)) THEN
+              aboveind(i,j,irm) = iv
+              newdist(i,j,iv) = varvalues(iv) + dTOTamount
+              irm = irm + 1
+            END IF
+          END DO
+
+          ! Final total amount of the new distribution
+          eamount = SUM(newdist(i,j,:))
+
+          IF (dbg) THEN
+            Nredis = COUNT(redistributed(i,j,:) /= -1)
+            PRINT *,TRIM(bspc), Nredis, 'Values removed from ranges ________'
+            DO ia=1, Nredis
+              ind = redistributed(i,j,ia)
+              PRINT *,TRIM(bspc) // '    ', ia, ' prev:', varvalues(ind), 'new:', newdist(i,j,ind)
+            END DO
+            PRINT *,TRIM(bspc), NaboveHquants, 'Values above ________'
+            DO ia=1, NaboveHquants
+              ind = aboveind(i,j,ia)
+              PRINT *,TRIM(bspc) // '    ', ia, 'prev:', varvalues(ind), 'new:', newdist(i,j,ind)
+            END DO
+          END IF
+
+          ! Check
+          IF (ABS(iamount - eamount) > 1.e-5) THEN
+            PRINT *, warnmsg
+            PRINT *,'  ' // TRIM(fname) // ': final quantity:', eamount, ' &', iamount, ' does not' // &
+              ' fit: abs(iamount - eamount) > 1.e-5:', abs(iamount-eamount), ' TOTamount:', TOTamount
+            STOP
+          END IF
+
+          IF (dbg) THEN
+            varvalues = newdist(i,j,:)
+            ndist = MINVAL(varvalues)
+            xdist = MAXVAL(varvalues)
+            meandist = SUM(varvalues)/(d3*one)
+
+            IF (ALLOCATED(newquants)) DEALLOCATE(newquants)
+            ALLOCATE(newquants(Nquants+1), STAT=ierr)
+            msg = "Problem allocating 'newquants'"
+            CALL ErrMsg(msg, fname, ierr)
+  
+            CALL quantilesR_K(d3, varvalues, Nquants+1, newquants)
+
+            PRINT *,TRIM(bspc) // '  new distribution values _______'
+            PRINT *,TRIM(bspc) // '  final total amount:', eamount
+            PRINT *,TRIM(bspc) // '  min:', ndist, 'max:', xdist, 'mean:', meandist
+            PRINT *,TRIM(bspc) // '  initial quantiles:'
+            CALL PrintQuantilesR_K(20, varvalues, Nquants+1, quants, Nstrings(spc,8))
+            PRINT *,TRIM(bspc) // '  new quantiles:'
+            CALL PrintQuantilesR_K(20, varvalues, Nquants+1, newquants, Nstrings(spc,8))
+            !STOP
+          END IF
+
+        ELSE
+
+          PRINT *,TRIM(warnmsg)
+          PRINT *,'  ' // TRIM(fname) // ': Any value has been removed !!'
+          Habovequant(i,j) = zero
+
+        END IF NonullRM
+
+      END IF doingnothing
+    END DO
+    PRINT *,'  ' // TRIM(fname) // ': done', i, ' of', d1
+  END DO
+
+  IF (ALLOCATED(varvalues)) DEALLOCATE(varvalues)
+  IF (ALLOCATED(rngind)) DEALLOCATE(rngind)
+  IF (ALLOCATED(rmind)) DEALLOCATE(rmind)
+  IF (ALLOCATED(quants)) DEALLOCATE(quants)
+  IF (ALLOCATED(newquants)) DEALLOCATE(newquants)
+
+  RETURN
+
+END SUBROUTINE randLowRemovePercent_HighEqui3D
+
+END MODULE module_DistriCorrection
Index: /trunk/tools/module_definitions.f90
===================================================================
--- /trunk/tools/module_definitions.f90	(revision 1609)
+++ /trunk/tools/module_definitions.f90	(revision 1609)
@@ -0,0 +1,16 @@
+MODULE module_definitions
+
+  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
+  REAL(r_k), PARAMETER                                   :: zeroRK = 0.
+  REAL(r_k), PARAMETER                                   :: oneRK = 1.
+  REAL(r_k), PARAMETER                                   :: twoRK = 2.
+  CHARACTER(len=100)                                     :: infomsg = 'INFORMATION -- information' // &
+    ' -- INFORMATION -- information'
+  CHARACTER(len=50)                                      :: warnmsg = 'WARNING -- warning -- ' //     &
+    'WARNING -- warning'
+  CHARACTER(len=50)                                      :: emsg = 'ERROR -- error -- ERROR -- error'
+  ! Fill value at 64 bits
+  REAL(r_k)                                              :: fillval64 = 1.e20
+  CHARACTER(len=50)                                      :: fname
+
+END MODULE module_definitions
Index: /trunk/tools/module_scientific.f90
===================================================================
--- /trunk/tools/module_scientific.f90	(revision 1609)
+++ /trunk/tools/module_scientific.f90	(revision 1609)
@@ -0,0 +1,250 @@
+MODULE module_scientific
+! Module of the scientific function/subroutines
+  
+!!!!!!! Functions & subroutines
+! FindMinimumR_K*: Function returns the location of the minimum in the section between Start and End.
+! PrintQuantilesR_K: Subroutine to print the quantiles of values REAL(r_k)
+! quantilesR_K: Subroutine to provide the quantiles of a given set of values of type real 'r_k'
+! rand_sample: Subroutine to randomly sample a range of indices
+! SortR_K*: Subroutine receives an array x() r_K and sorts it into ascending order.
+! SwapR_K*: Subroutine swaps the values of its two formal arguments.
+
+!!! *Functions/Subroutines to sort values adpated. The method used is usually referred to as "selection" method.
+! from: http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90
+
+  USE module_definitions
+  USE module_generic
+
+  CONTAINS
+
+SUBROUTINE rand_sample(Nvals, Nsample, sample)
+! Subroutine to randomly sample a range of indices
+
+  IMPLICIT NONE
+
+  INTEGER, INTENT(in)                                    :: Nvals, Nsample
+  INTEGER, DIMENSION(Nsample), INTENT(out)               :: sample
+
+! Local
+  INTEGER                                                :: i, ind, jmax
+  REAL, DIMENSION(Nsample)                               :: randv
+  CHARACTER(len=50)                                      :: fname
+  LOGICAL                                                :: found
+  LOGICAL, DIMENSION(Nvals)                              :: issampled
+  CHARACTER(len=256)                                     :: msg
+  CHARACTER(len=10)                                      :: IS1, IS2
+
+!!!!!!! Variables
+! Nvals: number of values
+! Nsamples: number of samples
+! sample: samnple 
+  fname = 'rand_sample'
+
+  IF (Nsample > Nvals) THEN
+    WRITE(IS1,'(I10)')Nvals
+    WRITE(IS2,'(I10)')Nsample
+    msg = 'Sampling of ' // TRIM(IS1) // ' is too big for ' // TRIM(IS1) // 'values'
+    CALL ErrMsg(msg, fname, -1)
+  END IF
+
+  ! Generation of random numbers always the same series during the whole program!
+  CALL RANDOM_NUMBER(randv)
+
+  ! Making sure that we do not repeat any value
+  issampled = .FALSE.
+
+  DO i=1, Nsample
+    ! Generation of the index from the random numbers
+    ind = MAX(INT(randv(i)*Nvals), 1)
+
+    IF (.NOT.issampled(ind)) THEN
+      sample(i) = ind
+      issampled(ind) = .TRUE.
+    ELSE
+      ! Looking around the given index
+      !PRINT *,' Index :', ind, ' already sampled!', issampled(ind)
+      found = .FALSE.
+      DO jmax=1, Nvals
+        ind = MIN(ind+jmax, Nvals)
+        IF (.NOT.issampled(ind)) THEN
+          sample(i) = ind
+          issampled(ind) = .TRUE.
+          found = .TRUE.
+          EXIT
+        END IF
+        ind = MAX(1, ind-jmax)
+        IF (.NOT.issampled(ind)) THEN
+          sample(i) = ind
+          issampled(ind) = .TRUE.
+          found = .TRUE.
+          EXIT
+        END IF
+      END DO
+      IF (.NOT.found) THEN
+        msg = 'sampling could not be finished due to absence of available value!!'
+        CALL ErrMsg(msg, fname, -1)
+      END IF
+    END IF
+
+  END DO
+
+  RETURN
+
+END SUBROUTINE rand_sample
+
+SUBROUTINE PrintQuantilesR_K(Nvals, vals, Nquants, qtvs, bspc)
+! Subroutine to print the quantiles of values REAL(r_k)
+
+  IMPLICIT NONE
+
+  INTEGER, INTENT(in)                                    :: Nvals, Nquants
+  REAL(r_k), DIMENSION(Nvals), INTENT(in)                :: vals
+  REAL(r_k), DIMENSION(Nquants), INTENT(in)              :: qtvs
+  CHARACTER(len=1000), OPTIONAL                          :: bspc
+
+! Local
+  INTEGER                                                :: iq
+  LOGICAL, DIMENSION(Nvals)                              :: search1, search2, search
+  CHARACTER(len=6)                                       :: RS1
+  CHARACTER(len=50)                                      :: fname
+  CHARACTER(len=1000)                                    :: bspcS
+
+!!!!!!! Variables
+! vals: series of values
+! qtvs: values of the quantiles
+! bspc: base quantity of spaces
+
+  fname = 'PrintQuantilesR_K'
+
+  IF (PRESENT(bspc)) THEN
+    bspcS = bspc
+  ELSE
+    bspcS = '      '
+  END IF
+
+  DO iq=1, Nquants-1
+
+    WHERE (vals >= qtvs(iq)) 
+      search1 = .TRUE.
+    ELSEWHERE
+      search1 = .FALSE.
+    END WHERE
+
+    WHERE (vals < qtvs(iq+1)) 
+      search2 = .TRUE.
+    ELSEWHERE
+      search2 = .FALSE.
+    END WHERE
+
+    WHERE (search1 .AND. search2) 
+      search = .TRUE.
+    ELSEWHERE
+      search = .FALSE.
+    END WHERE
+
+    WRITE(RS1, '(F6.2)')(iq)*100./(Nquants-1)
+    PRINT *, TRIM(bspcS) // '[',iq,']', TRIM(RS1) // ' %:', qtvs(iq), 'N:', COUNT(search)
+
+  END DO
+
+  RETURN
+
+END SUBROUTINE PrintQuantilesR_K
+
+   INTEGER FUNCTION FindMinimumR_K(x, dsize, Startv, Endv)
+! Function returns the location of the minimum in the section between Start and End.
+
+      IMPLICIT NONE
+
+      INTEGER, INTENT(in)                                :: dsize
+      REAL(r_k), DIMENSION(dsize), INTENT(in)            :: x
+      INTEGER, INTENT(in)                                :: Startv, Endv
+
+! Local
+      REAL(r_k)                                          :: Minimum
+      INTEGER                                            :: Location
+      INTEGER                                            :: i
+
+      Minimum  = x(Startv)                               ! assume the first is the min
+      Location = Startv                                  ! record its position
+      DO i = Startv+1, Endv                              ! start with next elements
+         IF (x(i) < Minimum) THEN                        !   if x(i) less than the min?
+            Minimum  = x(i)                              !      Yes, a new minimum found
+            Location = i                                 !      record its position
+         END IF
+      END DO
+
+      FindMinimumR_K = Location                          ! return the position
+
+   END FUNCTION  FindMinimumR_K
+
+   SUBROUTINE SwapR_K(a, b)
+! Subroutine swaps the values of its two formal arguments.
+
+      IMPLICIT  NONE
+
+      REAL(r_k), INTENT(INOUT)                           :: a, b
+! Local
+      REAL(r_k)                                          :: Temp
+
+      Temp = a
+      a    = b
+      b    = Temp
+
+   END SUBROUTINE  SwapR_K
+
+   SUBROUTINE  SortR_K(x, Nx)
+! Subroutine receives an array x() r_K and sorts it into ascending order.
+
+      IMPLICIT NONE
+
+      INTEGER, INTENT(IN)                                :: Nx
+      REAL(r_k), DIMENSION(Nx), INTENT(INOUT)            :: x
+
+! Local
+      INTEGER                                            :: i
+      INTEGER                                            :: Location
+
+      DO i = 1, Nx-1                                     ! except for the last
+         Location = FindMinimumR_K(x, Nx-i+1, i, Nx)     ! find min from this to last
+         CALL  SwapR_K(x(i), x(Location))                ! swap this and the minimum
+      END DO
+
+   END SUBROUTINE  SortR_K
+
+SUBROUTINE quantilesR_K(Nvals, vals, Nquants, quants)
+! Subroutine to provide the quantiles of a given set of values of type real 'r_k'
+
+  IMPLICIT NONE
+
+  INTEGER, INTENT(in)                                    :: Nvals, Nquants
+  REAL(r_k), DIMENSION(Nvals), INTENT(in)                :: vals
+  REAL(r_k), DIMENSION(Nquants), INTENT(out)             :: quants
+
+! Local
+  INTEGER                                                :: i
+  REAL(r_k)                                              :: minv, maxv
+  REAL(r_k), DIMENSION(Nvals)                            :: sortedvals
+
+!!!!!!! Variables
+! Nvals: number of values
+! Rk: kind of real of the values
+! vals: values
+! Nquants: number of quants
+! quants: values at which the quantile start
+
+  minv = MINVAL(vals)
+  maxv = MAXVAL(vals)
+
+  sortedvals = vals
+  ! Using from: http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90
+  CALL SortR_K(sortedvals, Nvals)
+
+  quants(1) = minv
+  DO i=2, Nquants
+    quants(i) = sortedvals(INT((i-1)*Nvals/Nquants))
+  END DO
+
+END SUBROUTINE quantilesR_K
+
+END MODULE module_scientific
