Changeset 2235 in lmdz_wrf

Timestamp:

Nov 16, 2018, 10:27:45 PM (6 years ago)

Author:

lfita

Message:

Working version of `temporal_stats' with 'percentiles'

Location:

trunk/tools

Files:

• module_generic.f90 (modified) (2 diffs)
• module_scientific.f90 (modified) (2 diffs)
• nc_var_tools.py (modified) (6 diffs)

trunk/tools/module_generic.f90

@@ -13 +13 @@
 ! Nvalues_2DArrayI: Number of different values of a 2D integer array
 ! mat2DPosition: Function to provide the i, j indices of a given value inside a 2D matrix
+! numberTimes: Function to provide the number of times that a given set of characters happen within a string
 ! RangeI: Function to provide a range of d1 values from 'iniv' to 'endv', of integer values in a vector
 ! RangeR: Function to provide a range of d1 values from 'iniv' to 'endv', of real values in a vector
@@ -225 +226 @@
 
   END FUNCTION Index2DArrayR_K
+
+  INTEGER FUNCTION numberTimes(String, charv)
+! Function to provide the number of times that a given set of characters happen within a string
+
+    IMPLICIT NONE
+
+    CHARACTER(LEN=*), INTENT(IN)                         :: String, charv
+
+! Local
+    INTEGER                                              :: i, Lstring, Lcharv
+
+    numberTimes = 0
+
+    Lstring = LEN_TRIM(String)
+    Lcharv = LEN_TRIM(charv)
+
+    DO i=1,Lstring - Lcharv
+      IF (String(i:i+Lcharv-1) == TRIM(charv)) numberTimes = numberTimes + 1
+    END DO
+
+  END FUNCTION numberTimes
+
 
   FUNCTION RangeI(d1, iniv, endv)

trunk/tools/module_scientific.f90

@@ -24 +24 @@
 ! paths_border: Subroutine to search the paths of a border field.
 ! path_properties: Subroutine to determine the properties of a path
+! percentiles_R_K*D: Subroutine to compute the percentiles of a *D R_K array along given set of axis
 ! polygon_properties: Subroutine to determine the properties of a polygon (as .TRUE. matrix)
 ! polygons: Subroutine to search the polygons of a border field. FORTRAN based. 1st = 1!
@@ -3682 +3683 @@
   END SUBROUTINE runmean_F1D
 
+  SUBROUTINE percentiles_R_K2D(values, axisS, Npercen, d1, d2, percentiles)
+  ! Subroutine to compute the percentiles of a 2D R_K array along given set of axis
+
+    IMPLICIT NONE
+  
+    INTEGER, INTENT(in)                                    :: d1, d2, Npercen
+    REAL(r_k), DIMENSION(d1,d2), INTENT(in)                :: values
+    CHARACTER(LEN=*), INTENT(in)                           :: axisS
+    REAL(r_k), DIMENSION(d1, d2, Npercen), INTENT(out)     :: percentiles
+
+    ! Local
+    INTEGER                                                :: i
+    INTEGER                                                :: Lstring, LaxisS, iichar
+    CHARACTER(LEN=1000)                                    :: splitaxis
+    INTEGER, DIMENSION(1)                                  :: axis1
+    CHARACTER(LEN=200), DIMENSION(2)                       :: axis2S
+    INTEGER, DIMENSION(2)                                  :: axis2
+    CHARACTER(LEN=1)                                       :: Naxs
+
+!!!!!!! Variables
+! d1,d2: length of the 2D dimensions
+! values: values to use to compute the percentiles
+! axisS: ':' separated list of axis to use to compute the percentiles ('all' for all axes)
+! Npercen: number of percentiles
+! percentiles: percentiles of the daata
+
+    fname = 'percentiles_R_K2D'
+
+    LaxisS = LEN_TRIM(axisS)
+    iichar = numberTimes(axisS(1:LaxisS), ':')
+
+    splitaxis = ''
+    splitaxis(1:LaxisS) = axisS(1:LaxisS)
+    percentiles = 0.
+
+    IF (iichar == 0) THEN
+      READ(axisS,'(I1)')axis1(1)
+    ELSE IF (iichar == 1) THEN
+      CALL split(splitaxis, ':', 2, axis2S)
+    ELSE
+      WRITE(Naxs,'(A1)')iichar
+      msg = "' rank 2 values can not compute percentiles using " // Naxs // "' number of axis !!"
+      CALL ErrMsg(msg, fname, -1)
+    END IF
+
+    IF (TRIM(axisS) == 'all') iichar = 2
+  
+    IF (iichar == 0) THEN
+      ! Might be a better way, but today I can't think it !!
+      IF (axis1(1) == 1) THEN
+        DO i=1, d2
+          CALL quantilesR_K(d1, values(:,i), Npercen, percentiles(1,i,:))
+        END DO
+      ELSE IF (axis1(1) == 2) THEN
+        DO i=1, d1
+          CALL quantilesR_K(d2, values(i,:), Npercen, percentiles(i,1,:))
+        END DO
+      ELSE
+        WRITE(Naxs,'(A1)')axis1(1)
+        msg = "' rank 2 values can not compute percentiles using axis " // Naxs // "' !!"
+        CALL ErrMsg(msg, fname, -1)
+      END IF
+    ELSE
+      CALL quantilesR_K(d1*d2, RESHAPE(values, (/d1*d2/)), Npercen, percentiles(1,1,:))
+    END IF
+
+  END SUBROUTINE percentiles_R_K2D
+
+  SUBROUTINE percentiles_R_K3D(values, axisS, Npercen, d1, d2, d3, percentiles)
+  ! Subroutine to compute the percentiles of a 3D R_K array along given set of axis
+
+    IMPLICIT NONE
+  
+    INTEGER, INTENT(in)                                    :: d1, d2, d3, Npercen
+    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)             :: values
+    CHARACTER(LEN=*), INTENT(in)                           :: axisS
+    REAL(r_k), DIMENSION(d1, d2, d3, Npercen), INTENT(out) :: percentiles
+
+    ! Local
+    INTEGER                                                :: i, j
+    INTEGER                                                :: Lstring, LaxisS, iichar
+    CHARACTER(LEN=1000)                                    :: splitaxis
+    INTEGER, DIMENSION(1)                                  :: axis1
+    CHARACTER(LEN=200), DIMENSION(2)                       :: axis2S
+    INTEGER, DIMENSION(2)                                  :: axis2
+    CHARACTER(LEN=200), DIMENSION(3)                       :: axis3S
+    INTEGER, DIMENSION(3)                                  :: axis3
+    CHARACTER(LEN=1)                                       :: Naxs1, Naxs2
+
+!!!!!!! Variables
+! d1,d2: length of the 2D dimensions
+! values: values to use to compute the percentiles
+! axisS: ':' separated list of axis to use to compute the percentiles ('all' for all axes)
+! Npercen: number of percentiles
+! percentiles: percentiles of the daata
+
+    fname = 'percentiles_R_K3D'
+
+    LaxisS = LEN_TRIM(axisS)
+    iichar = numberTimes(axisS(1:LaxisS), ':')
+
+    splitaxis = ''
+    splitaxis(1:LaxisS) = axisS(1:LaxisS)
+
+    percentiles = 0.
+
+    IF (iichar == 0) THEN
+      READ(axisS,'(I1)')axis1(1)
+    ELSE IF (iichar == 1) THEN
+      CALL split(splitaxis, ':', 2, axis2S)
+      DO i=1,2
+        READ(axis2S(i), '(I1)')axis2(i)
+      END DO
+    ELSE IF (iichar == 2) THEN
+      CALL split(splitaxis, ':', 3, axis3S)
+    ELSE
+      READ(Naxs1,'(A1)')iichar
+      msg = "' rank 3 values can not compute percentiles using " // Naxs1 // "' number of axis !!"
+      CALL ErrMsg(msg, fname, -1)
+    END IF
+
+    IF (TRIM(axisS) == 'all') iichar = 3
+  
+    IF (iichar == 0) THEN
+      ! Might be a better way, but today I can't think it !!
+      IF (axis1(1) == 1) THEN
+        DO i=1, d2
+          DO j=1, d3
+            CALL quantilesR_K(d1, values(:,i,j), Npercen, percentiles(1,i,j,:))
+          END DO
+        END DO
+      ELSE IF (axis1(1) == 2) THEN
+        DO i=1, d1
+          DO j=1, d3
+            CALL quantilesR_K(d2, values(i,:,j), Npercen, percentiles(i,1,j,:))
+          END DO
+        END DO
+      ELSE IF (axis1(1) == 3) THEN
+        DO i=1, d1
+          DO j=1, d2
+            CALL quantilesR_K(d3, values(i,j,:), Npercen, percentiles(i,j,1,:))
+          END DO
+        END DO
+      ELSE
+        WRITE(Naxs1,'(A1)')axis1(1)
+        msg = "' rank 3 values can not compute percentiles using axis " // Naxs1 // "' !!"
+        CALL ErrMsg(msg, fname, -1)
+      END IF
+    ELSE IF (iichar == 1) THEN
+      ! Might be a better way, but today I can't think it !!
+      IF (axis2(1) == 1 .AND. axis2(2) == 2) THEN
+        DO i=1, d3
+          CALL quantilesR_K(d1*d2, RESHAPE(values(:,:,i), (/d1*d2/)), Npercen, percentiles(1,1,i,:))
+        END DO
+      ELSE IF (axis2(1) == 1 .AND. axis2(2) == 3) THEN
+        DO i=1, d2
+          CALL quantilesR_K(d1*d3, RESHAPE(values(:,i,:), (/d1*d3/)), Npercen, percentiles(1,i,1,:))
+        END DO
+      ELSE IF (axis2(1) == 2 .AND. axis2(2) == 3) THEN
+        DO i=1, d1
+          CALL quantilesR_K(d2*d3, RESHAPE(values(i,:,:), (/d2*d3/)), Npercen, percentiles(i,1,1,:))
+        END DO
+      ELSE
+        WRITE(Naxs1,'(A1)')axis2(1)
+        WRITE(Naxs2,'(A1)')axis2(2)
+        msg="' rank 3 values can not compute percentiles using axis "//Naxs1// ', ' // Naxs2 // "' !!"
+        CALL ErrMsg(msg, fname, -1)
+      END IF
+    ELSE
+      CALL quantilesR_K(d1*d2*d3, RESHAPE(values, (/d1*d2*d3/)), Npercen, percentiles(1,1,1,:))
+    END IF
+
+  END SUBROUTINE percentiles_R_K3D
+
+  SUBROUTINE percentiles_R_K4D(values, axisS, Npercen, d1, d2, d3, d4, percentiles)
+  ! Subroutine to compute the percentiles of a 4D R_K array along given set of axis
+
+    IMPLICIT NONE
+  
+    INTEGER, INTENT(in)                                    :: d1, d2, d3, d4, Npercen
+    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)          :: values
+    CHARACTER(LEN=*), INTENT(in)                           :: axisS
+    REAL(r_k), DIMENSION(d1,d2,d3,d4,Npercen), INTENT(out) :: percentiles
+
+    ! Local
+    INTEGER                                                :: i, j, k
+    INTEGER                                                :: Lstring, LaxisS, iichar
+    CHARACTER(LEN=1000)                                    :: splitaxis
+    INTEGER, DIMENSION(1)                                  :: axis1
+    CHARACTER(LEN=200), DIMENSION(2)                       :: axis2S
+    INTEGER, DIMENSION(2)                                  :: axis2
+    CHARACTER(LEN=200), DIMENSION(3)                       :: axis3S
+    INTEGER, DIMENSION(3)                                  :: axis3
+    CHARACTER(LEN=200), DIMENSION(4)                       :: axis4S
+    CHARACTER(LEN=1)                                       :: Naxs1, Naxs2, Naxs3
+
+!!!!!!! Variables
+! d1,d2: length of the 2D dimensions
+! values: values to use to compute the percentiles
+! axisS: ':' separated list of axis to use to compute the percentiles ('all' for all axes)
+! Npercen: number of percentiles
+! percentiles: percentiles of the daata
+
+    fname = 'percentiles_R_K3D'
+
+    LaxisS = LEN_TRIM(axisS)
+    iichar = numberTimes(axisS(1:LaxisS), ':')
+
+    splitaxis = ''
+    splitaxis(1:LaxisS) = axisS(1:LaxisS)
+
+    percentiles = 0.
+
+    PRINT *,'iichar:', iichar, axisS(1:LaxisS)
+
+    IF (iichar == 0) THEN
+      READ(axisS,'(I1)')axis1(1)
+    ELSE IF (iichar == 1) THEN
+      CALL split(splitaxis, ':', 2, axis2S)
+      DO i=1,2
+        READ(axis2S(i), '(I1)')axis2(i)
+      END DO
+    ELSE IF (iichar == 2) THEN
+      CALL split(splitaxis, ':', 3, axis3S)
+      DO i=1,3
+        READ(axis3S(i), '(I1)')axis3(i)
+      END DO
+    ELSE IF (iichar == 3) THEN
+      CALL split(splitaxis, ':', 4, axis4S)
+    ELSE
+      READ(Naxs1,'(A1)')iichar
+      msg = "' rank 4 values can not compute percentiles using " // Naxs1 // "' number of axis !!"
+      CALL ErrMsg(msg, fname, -1)
+    END IF
+
+    IF (TRIM(axisS) == 'all') iichar = 4
+  
+    IF (iichar == 0) THEN
+      ! Might be a better way, but today I can't think it !!
+      IF (axis1(1) == 1) THEN
+        DO i=1, d2
+          DO j=1, d3
+            DO k=1, d4
+              CALL quantilesR_K(d1, values(:,i,j,k), Npercen, percentiles(1,i,j,k,:))
+            END DO
+          END DO
+        END DO
+      ELSE IF (axis1(1) == 2) THEN
+        DO i=1, d1
+          DO j=1, d3
+            DO k=1, d4
+              CALL quantilesR_K(d2, values(i,:,j,k), Npercen, percentiles(i,1,j,k,:))
+            END DO
+          END DO
+        END DO
+      ELSE IF (axis1(1) == 3) THEN
+        DO i=1, d1
+          DO j=1, d2
+            DO k=1, d4
+              CALL quantilesR_K(d3, values(i,j,:,k), Npercen, percentiles(i,j,1,k,:))
+            END DO
+          END DO
+        END DO
+      ELSE IF (axis1(1) == 4) THEN
+        DO i=1, d1
+          DO j=1, d2
+            DO k=1, d3
+              CALL quantilesR_K(d4, values(i,j,k,:), Npercen, percentiles(i,j,k,1,:))
+            END DO
+          END DO
+        END DO
+      ELSE
+        WRITE(Naxs1,'(A1)')axis1(1)
+        msg = "' rank 3 values can not compute percentiles using axis " // Naxs1 // "' !!"
+        CALL ErrMsg(msg, fname, -1)
+      END IF
+    ELSE IF (iichar == 1) THEN
+      ! Might be a better way, but today I can't think it !!
+      IF (axis2(1) == 1 .AND. axis2(2) == 2) THEN
+        DO i=1, d3
+          DO j=1, d4
+            CALL quantilesR_K(d1*d2, RESHAPE(values(:,:,i,j), (/d1*d2/)), Npercen,                    &
+              percentiles(1,1,i,j,:))
+          END DO
+        END DO
+      ELSE IF (axis2(1) == 1 .AND. axis2(2) == 3) THEN
+        DO i=1, d2
+          DO j=1, d4
+            CALL quantilesR_K(d1*d3, RESHAPE(values(:,i,:,j), (/d1*d3/)), Npercen,                    &
+              percentiles(1,i,1,j,:))
+          END DO
+        END DO
+      ELSE IF (axis2(1) == 1 .AND. axis2(2) == 4) THEN
+        DO i=1, d2
+          DO j=1, d3
+            CALL quantilesR_K(d1*d4, RESHAPE(values(:,i,j,:), (/d1*d4/)), Npercen,                    &
+              percentiles(1,i,j,1,:))
+          END DO
+        END DO
+      ELSE IF (axis2(1) == 2 .AND. axis2(2) == 3) THEN
+        DO i=1, d1
+          DO j=1, d4
+            CALL quantilesR_K(d2*d3, RESHAPE(values(i,:,:,j), (/d2*d3/)), Npercen,                    &
+              percentiles(i,1,1,j,:))
+          END DO
+        END DO
+      ELSE IF (axis2(1) == 2 .AND. axis2(2) == 4) THEN
+        DO i=1, d1
+          DO j=1, d3
+            CALL quantilesR_K(d2*d4, RESHAPE(values(i,:,j,:), (/d2*d4/)), Npercen,                    &
+              percentiles(i,1,j,1,:))
+          END DO
+        END DO
+      ELSE IF (axis2(1) == 3 .AND. axis2(2) == 4) THEN
+        DO i=1, d1
+          DO j=1, d2
+            CALL quantilesR_K(d3*d4, RESHAPE(values(i,j,:,:), (/d3*d4/)), Npercen,                    &
+              percentiles(i,j,1,1,:))
+          END DO
+        END DO
+      ELSE
+        WRITE(Naxs1,'(A1)')axis2(1)
+        WRITE(Naxs2,'(A1)')axis2(2)
+        msg="' rank 4 values can not compute percentiles using axis "//Naxs1// ', ' // Naxs2 // "' !!"
+        CALL ErrMsg(msg, fname, -1)
+      END IF
+    ELSE IF (iichar == 2) THEN
+      IF (axis2(1) == 1 .AND. axis2(2) == 2 .AND. axis3(3) == 3) THEN
+        DO i=1, d4
+          CALL quantilesR_K(d1*d2*d3, RESHAPE(values(:,:,:,i), (/d1*d2*d3/)), Npercen,                &
+            percentiles(1,1,1,i,:))
+        END DO
+      ELSE IF (axis2(1) == 1 .AND. axis2(2) == 2 .AND. axis3(3) == 4) THEN
+        DO i=1, d3
+          CALL quantilesR_K(d1*d2*d4, RESHAPE(values(:,:,i,:), (/d1*d2*d4/)), Npercen,                &
+            percentiles(1,1,i,1,:))
+        END DO
+      ELSE IF (axis2(1) == 1 .AND. axis2(2) == 3 .AND. axis3(3) == 4) THEN
+        DO i=1, d2
+          CALL quantilesR_K(d1*d3*d4, RESHAPE(values(:,i,:,:), (/d1*d3*d4/)), Npercen,                &
+            percentiles(1,i,1,1,:))
+        END DO
+      ELSE IF (axis2(1) == 2 .AND. axis2(2) == 3 .AND. axis3(3) == 4) THEN
+        DO i=1, d1
+          CALL quantilesR_K(d2*d3*d4, RESHAPE(values(i,:,:,:), (/d2*d3*d4/)), Npercen,                &
+            percentiles(i,1,1,1,:))
+        END DO
+      ELSE
+        WRITE(Naxs1,'(A1)')axis3(1)
+        WRITE(Naxs2,'(A1)')axis3(2)
+        WRITE(Naxs3,'(A1)')axis3(2)
+        msg="' rank 4 values can not compute percentiles using axis "// Naxs1 // ', ' // Naxs2 //     &
+          ', ' // Naxs3 //"' !!"
+        CALL ErrMsg(msg, fname, -1)
+      END IF
+    ELSE
+      CALL quantilesR_K(d1*d2*d3*d4, RESHAPE(values, (/d1*d2*d3*d4/)), Npercen, percentiles(1,1,1,1,:))
+    END IF
+
+  END SUBROUTINE percentiles_R_K4D
+
 END MODULE module_scientific

trunk/tools/nc_var_tools.py

@@ -17 +17 @@
 import subprocess as sub
 import module_ForDef as fdef
+import module_ForSci as fsci
 
 ####### Contents:
@@ -25431 +25432 @@
         timestatsv.append(amount)
         statn = timest.split('@')[iamount+1]
+        if not gen.searchInlist(Lstatn.keys(), statn):
+            print errormsg
+            print '  ' + fname + ": statistics '" + statn + "' not ready !!"
+            print '    available ones:', Lstatn.keys()
+            quit(-1)
+
         timestatsv.append(statn)
-        if statn == 'percentiles': 
+        if statn == 'percentiles':
+            if len(timest.split('@')) != iamount+3:
+                print errormsg
+                print '  ' + fname + ': no number of percentiles provided!!'
+                print '    provided values:', timest
+                quit(-1)
             Npercen = int(timest.split('@')[iamount+2])
             timestatsv.append(Npercen)
@@ -25464 +25476 @@
         amount = timestatv[iamount]
 
-        timeslice, Ntimeslice = gen.time_slices(timevc, timeu, timec, period, amount)
+        time_desc = gen.temporal_desc(timevc, timeu, timec)
+        timeslice, Ntimeslice= gen.time_slices(timevc, timeu, timec, period, amount, \
+          time_desc)
         statslices[iTst] = [Ntimeslice, timeslice]
  
@@ -25600 +25614 @@
                 newvarn = varn+str(amount)+periodn+statn
 
-            newvarv = onewnc.createVariable(newvarn, 'f', tuple(dimntstat),         \
-              fill_value=gen.fillValueR)
+            if statn != 'percentiles':
+                newvarv = onewnc.createVariable(newvarn, 'f', tuple(dimntstat),     \
+                  fill_value=gen.fillValueF)
+            else:
+                if not gen.searchInlist(onewnc.dimensions, 'percentile'):
+                    print '  ' + fname + ": adding '" + statn + "' dimension !!"
+                    newdim = onewnc.createDimension('percentile', Npercen+1)
+                
+                    newvar= onewnc.createVariable('percentile','f', ('percentile'))
+                    newvar[:] = np.arange(0.,100.+100./Npercen,100./Npercen)
+                    basicvardef(newvar, 'percentile', 'Percentile', '%')
+
+                dimntstat = ['percentile'] + dimntstat
+                newvarv = onewnc.createVariable(newvarn, 'f', tuple(dimntstat),     \
+                  fill_value=gen.fillValueF)
+                
             varunits = get_varunits(ovar)
             if statn == 'std': varunits = '(' + varunits + ')**2)'
@@ -25624 +25652 @@
                 newvarN[islc] = timeslcv[1]-timeslcv[0]+1
                 
-                timeslc = []
-                shapetimeslc = []
-                tstatslc = []
-                for dn in vardims:
-                    if dn == timedimn:
-                        timeslc.append(slice(timeslcv[0],timeslcv[1],timeslcv[2]))
-                        shapetimeslc.append(timeslcv[1]-timeslcv[0]+1)
-                        tstatslc.append(islc)
-                    else:
-                        Ldim = len(onc.dimensions[dn])
-                        timeslc.append(slice(0,Ldim))
-                        shapetimeslc.append(Ldim)
-                        tstatslc.append(slice(0,Ldim))
-
-                tvals = ovar[tuple(timeslc)]
+                if periodn[0:3] != 'agg':
+                    timeslc = []
+                    shapetimeslc = []
+                    tstatslc = []
+                    for dn in vardims:
+                        if dn == timedimn:
+                            timeslc.append(slice(timeslcv[0],timeslcv[1],timeslcv[2]))
+                            shapetimeslc.append(timeslcv[1]-timeslcv[0]+1)
+                            tstatslc.append(islc)
+                        else:
+                            Ldim = len(onc.dimensions[dn])
+                            timeslc.append(slice(0,Ldim))
+                            shapetimeslc.append(Ldim)
+                            tstatslc.append(slice(0,Ldim))
+                    tvals = ovar[tuple(timeslc)]
+                else:
+                    newshape = list(vals.shape)
+                    aggslc = tslc[islc]
+                    NNtslc = len(aggslc)
+                    newshape[itdim] = NNtslc
+                    tvals = np.zeros(tuple(newshape), dtype=vals.dtype)
+                    origshape = list(vals.shape)
+                    iishape = list(vals.shape)
+                    for iislc in range(NNtslc):
+                        origshape[itdim] = aggslc[iislc]
+                        iishape[itdim] = iislc
+                        tvals[tuple(iishape)] = vals[tuple(origshape)]
+
                 if statn == 'min':
                     newvarv[tuple(tstatslc)] = np.min(tvals, axis=itdim)
@@ -25649 +25690 @@
                 elif statn == 'std':
                     newvarv[tuple(tstatslc)] = np.std(tvals, axis=itdim)
+                elif statn == 'percentiles':
+                    axesS = str(len(tvals.shape)-itdim)
+                    if tvals.shape[itdim] < Npercen:
+                        tstatslc = [slice(0,Npercen+1)] + tstatslc
+                        print warnmsg
+                        print '  ' + fname + ': Not enough temporal data:',          \
+                          tvals.shape[itdim], ' to copmute:',Npercen, 'percentiles !!'
+                        dimspercen = [Npercen+1]
+                        for i in range(len(tvals.shape)):
+                            if i != itdim:
+                                dimspercen.append(tvals.shape[i])
+                        percens = np.ones(tuple(dimspercen), dtype=np.float)*gen.fillValueF
+                        newvarv[tuple(tstatslc)] = percens[:]
+                    else:
+                        tvalst = tvals.transpose()
+                        if len(tvals.shape) == 4:
+                            percenst = fsci.module_scientific.percentiles_r_k4d(tvalst,  \
+                              axesS, npercen=Npercen+1, d1=tvals.shape[3],               \
+                              d2=tvals.shape[2], d3=tvals.shape[1], d4=tvals.shape[0])
+                        elif len(tvals.shape) == 3:
+                            percenst = fsci.module_scientific.percentiles_r_k3d(tvalst,  \
+                              axesS, npercen=Npercen+1, d1=tvals.shape[2],               \
+                              d2=tvals.shape[1], d3=tvals.shape[0])
+                        elif len(tvals.shape) == 2:
+                            percenst = fsci.module_scientific.percentiles_r_k2d(tvalst,  \
+                              axesS, npercen=Npercen+1, d1=tvals.shape[1],               \
+                              d2=tvals.shape[0])
+                        elif len(tvals.shape) == 1:
+                            quants = Quantiles(tvals, Npercen)
+                            percenst = quants.quantilesv
+                        else:
+                            print errormsg
+                            print '  ' +fname+ ': rank of temporal data:', tvals.shape,  \
+                              'not ready !!'
+                            quit(-1)
+                        percens = percenst.transpose()
+
+                        shapev = list(percens.shape)
+                        for idd in range(len(shapev)):
+                            if idd != itdim+1:
+                                shapev[idd] = slice(0,shapev[idd])
+                            else:
+                                shapev[itdim+1] = 0
+
+                        tstatslc = [slice(0,Npercen+1)] + tstatslc
+                        newvarv[tuple(tstatslc)] = percens[tuple(shapev)]
                 else:
                     print errormsg