Changeset 2304 in lmdz_wrf for trunk

Timestamp:

Jan 30, 2019, 7:41:32 PM (6 years ago)

Author:

lfita

Message:

Almost working !?

Location:

trunk/tools

Files:

• module_scientific.f90 (modified) (11 diffs)
• nc_var_tools.py (modified) (7 diffs)

trunk/tools/module_scientific.f90

@@ -44 +44 @@
 ! multi_index_mat3DRK: Subroutine to provide the indices of the different locations of a value inside a 3D RK matrix
 ! multi_index_mat4DRK: Subroutine to provide the indices of the different locations of a value inside a 4D RK matrix
+! multi_spaceweightstats_in3DRK3_slc3v3: Subroutine to compute an spatial statistics value from a 3D RK matrix using 3rd dimension as 
+!   running one into a matrix of 3-variables slices of rank 3 using spatial weights
+! multi_spaceweightstats_in3DRK3_slc3v4: Subroutine to compute an spatial statistics value from a 3D RK matrix using 3rd dimension as 
+!   running one into a matrix of 3-variables slices of rank 4 using spatial weights
 ! NcountR: Subroutine to count real values
 ! paths_border: Subroutine to search the paths of a border field.
@@ -5634 +5638 @@
   END SUBROUTINE spaceweightstats
 
-  SUBROUTINE multi_spacewightstats_in3DRK3_slc3v3(varin, Ngridsin, gridsin, percentages, varout,      &
+  SUBROUTINE multi_spaceweightstats_in3DRK3_slc3v3(varin, Ngridsin, gridsin, percentages, varout,     &
     di1, di2, di3, ds1, ds2, ds3, maxNgridsin)
   ! Subroutine to compute an spatial statistics value from a 3D RK matrix using 3rd dimension as 
@@ -5653 +5657 @@
 ! Local
     INTEGER                                              :: i1, i2, i3, s1, s2, s3, iv
+    INTEGER                                              :: ii3, ss1, ss2, ss3
     INTEGER                                              :: Ncounts, Nin
     CHARACTER(len=3)                                     :: val1S, val2S
     CHARACTER(len=30)                                    :: val3S
     REAL(r_k)                                            :: minv, maxv, meanv, mean2v, stdv, medv
-    INTEGER, DIMENSION(:), ALLOCATABLE                   :: pin
+    REAL(r_k), DIMENSION(:), ALLOCATABLE                 :: pin
     INTEGER, DIMENSION(:,:), ALLOCATABLE                 :: gin
     REAL(r_k), DIMENSION(:), ALLOCATABLE                 :: svin
@@ -5681 +5686 @@
 ! varout: output statistical variable
 
-    fname = 'multi_spacewightstats_in3DRK3_slc3v3'
+    fname = 'multi_spaceweightstats_in3DRK3_slc3v3'
+
+    varout = fillval64
+
+    ss1 = 8 + 1
+    ss2 = 5 + 1
+    ss3 = 3 + 1
+    ii3 = 1 + 1
 
     ! Let's be efficient?
@@ -5700 +5712 @@
           pin = percentages(s1,s2,s3,1:Nin)
 
-          minv = fillVal64
-          maxv = -fillVal64
-          meanv = zeroRK
-          mean2v = zeroRK
-          stdv = zeroRK
-
           ! Getting the values
           DO iv=1, Nin
@@ -5713 +5719 @@
           END DO
 
+          IF (s1 == ss1 .AND. s2 == ss2 .AND. s3 == ss3) THEN
+            PRINT *, ' Lluis FOR Nin:', Nin
+            PRINT *, ' Lluis FOR values:', vin(:,ii3)
+            DO iv=1, Nin
+              PRINT *,' c:', iv, ':', gin(iv,:), ' p:', pin(iv)
+            END DO
+          END IF
+
           ! Computing along d3
           DO i3=1, di3
+            minv = fillVal64
+            maxv = -fillVal64
+            meanv = zeroRK
+            mean2v = zeroRK
+            stdv = zeroRK
             minv = MINVAL(vin(:,i3))
             maxv = MAXVAL(vin(:,i3))
-            meanv = SUM(vin(:,i3)*pin)/Nin
-            mean2v = SUM(vin(:,i3)**2*pin)/Nin
-            stdv = SQRT(mean2v - meanv**2)
+            meanv = SUM(vin(:,i3)*pin)
+            mean2v = SUM(vin(:,i3)**2*pin)
+            DO iv=1,Nin
+              stdv = stdv + ( (meanv - vin(iv,i3))*pin(iv) )**2
+            END DO
+            stdv = SQRT(stdv)
             svin = vin(:,i3)
             CALL SortR_K(svin, Nin)
@@ -5742 +5764 @@
     RETURN
 
-  END SUBROUTINE multi_spacewightstats_in3DRK3_slc3v3
-
-  SUBROUTINE multi_spacewightstats_in3DRK3_slc3v4(varin, Ngridsin, gridsin, percentages, varout,      &
+  END SUBROUTINE multi_spaceweightstats_in3DRK3_slc3v3
+
+  SUBROUTINE multi_spaceweightstats_in3DRK3_slc3v4(varin, Ngridsin, gridsin, percentages, varout,     &
     di1, di2, di3, ds1, ds2, ds3, ds4, maxNgridsin)
   ! Subroutine to compute an spatial statistics value from a 3D RK matrix using 3rd dimension as 
@@ -5768 +5790 @@
     CHARACTER(len=30)                                    :: val3S
     REAL(r_k)                                            :: minv, maxv, meanv, mean2v, stdv, medv
-    INTEGER, DIMENSION(:), ALLOCATABLE                   :: pin
+    REAL(r_k), DIMENSION(:), ALLOCATABLE                 :: pin
     INTEGER, DIMENSION(:,:), ALLOCATABLE                 :: gin
     REAL(r_k), DIMENSION(:), ALLOCATABLE                 :: svin
@@ -5792 +5814 @@
 ! varout: output statistical variable
 
-    fname = 'multi_spacewightstats_in3DRK3_slc3v4'
+    fname = 'multi_spaceweightstats_in3DRK3_slc3v4'
+
+    varout = fillval64
 
     ! Let's be efficient?
@@ -5827 +5851 @@
             ! Computing along d3
             DO i3=1, di3
+              minv = fillVal64
+              maxv = -fillVal64
+              meanv = zeroRK
+              mean2v = zeroRK
+              stdv = zeroRK
+
               minv = MINVAL(vin(:,i3))
               maxv = MAXVAL(vin(:,i3))
-              meanv = SUM(vin(:,i3)*pin)/Nin
-              mean2v = SUM(vin(:,i3)**2*pin)/Nin
-              stdv = SQRT(mean2v - meanv**2)
+              meanv = SUM(vin(:,i3)*pin)
+              mean2v = SUM(vin(:,i3)**2*pin)
+              DO iv=1,Nin
+                stdv = stdv + ( (meanv - vin(iv,i3))*pin(iv) )**2
+              END DO
+              stdv = SQRT(stdv)
               svin = vin(:,i3)
               CALL SortR_K(svin, Nin)
@@ -5855 +5888 @@
     RETURN
 
-  END SUBROUTINE multi_spacewightstats_in3DRK3_slc3v4
+  END SUBROUTINE multi_spaceweightstats_in3DRK3_slc3v4
 
   SUBROUTINE multi_index_mat2DI(d1, d2, d12, mat, value, Nindices, indices)

trunk/tools/nc_var_tools.py

@@ -28742 +28742 @@
                 for iB in range(dxB):
                     Nin = sN[jA,iA,jB,iB]
-                    sgpa[0:Nin,jA,iA,jB,iB] = anewvar[0:Nin,jA,iA,jB,iB]/            \
-                      aanewvar[jA,iA,jB,iB]
-                    panewvar[0:Nin,jA,iA,jB,iB] = sgpa[0:Nin,jA,iA,jB,iB]
+                    if Nin > 0 and aanewvar[jA,iA,jB,iB] != 0.:
+                        sgpa[0:Nin,jA,iA,jB,iB] = anewvar[0:Nin,jA,iA,jB,iB]/        \
+                          aanewvar[jA,iA,jB,iB]
+                        panewvar[0:Nin,jA,iA,jB,iB] = sgpa[0:Nin,jA,iA,jB,iB]
+                        # Let's check
+                        psum = np.sum(sgpa[0:Nin,jA,iA,jB,iB])
+                        if psum != 1.:
+                            print warnmsg
+                            print '    ' + fname + 'at slide:', jA,iA,jB,iB,         \
+                              'percentages do NOT add one !!'
+                            print '    Ngrids:', Nin, 'slice area:',                 \
+                              aanewvar[jA,iA,jB,iB], 'total sum:', psum
+                            #for iv in range(Nin):
+                            #    print iv, sgpa[iv,jA,iA,jB,iB]
     onewnc.sync()
                      
+    # Masking...
+    sgpa = ma.masked_equal(sgpa,gen.fillValueF)
+
     # removing monotones. Fortran only accepts rank <= 7 arrays !!
     sN, Ndims = gen.remove_monotones(sN, oNslice.dimensions)
@@ -28759 +28773 @@
 
     sNt = sN.transpose()
-    sint = sin.transpose()
+    #sint = sin.transpose()
+    sint0 = sin.transpose()
+    ## We are in Fortran thus (dx,dy) and 1st:1 !!!!
+    sint = np.ones((sint0.shape), dtype=np.float)*gen.fillValueF
+    sint[...,0] = sint0[...,1]
+    sint[...,1] = sint0[...,0]
     spt = sgpa.transpose()
+
+    # Checking
+    i = 3
+    j = 5
+    k = 8
+    chijk = [i,j,k]
+    it = 1
+    slicesvals = []
+    print 'Lluis Ndims:', Ndims
+    idd = 0
+    for dnn in Ndims:
+        scvv = slicesinf[dnn.split('_')[1]]
+        sccvvv = scvv[3]
+        print 'sccvvv:', sccvvv
+        print 'Lluis dnn:', dnn, sccvvv[chijk[idd]]
+        slicesvals.append(sccvvv[chijk[idd]])
+        idd = idd + 1
 
     for varn in varns:
@@ -28776 +28812 @@
         for dnn in vdimns:
             vdimn = dimvars[dnn]
-            if vdimn == 'WRFtime': 
+            if vdimn == 'WRFtime':
                 vdimn = 'time'
-                odimvar = onc.variables['Times']
-                timewrfv = odimvar[:]
-                tvals, urefvals = compute_WRFtime(timewrfv, refdate='19491201000000',\
-                  tunitsval='minutes')
-                tvals = np.array(tvals)
-                dimtwrf = len(tvals)
+                if not onewnc.variables.has_key('time'): 
+                    odimvar = onc.variables['Times']
+                    timewrfv = odimvar[:]
+                    tvals, urefvals = compute_WRFtime(timewrfv,                      \
+                      refdate='19491201000000', tunitsval='minutes')
+                    tvals = np.array(tvals)
+                    dimtwrf = len(tvals)
+                    if not gen.searchInlist(onewnc.dimensions,vdimn):
+                        newdim = onewnc.createDimension(vdimn, None)
+                    if not gen.searchInlist(onewnc.variables,vdimn):
+                        newvar = onewnc.createVariable(vdimn, 'f', ('time'))
+                        newvar[:] = tvals
+                        basicvardef(newvar, 'time', 'Time', urefvals)
+                        newvar.setncattr('calendar','gregorian')
                 if gen.searchInlist(slicestatsdim, dnn):
                     vardims.append(vdimn)
                     varshape.append(dimtwrf)
-                if not gen.searchInlist(onewnc.dimensions,vdimn):
-                    newdim = onewnc.createDimension(vdimn, None)
-                if not gen.searchInlist(onewnc.variables,vdimn):
-                    newvar = onewnc.createVariable(vdimn, 'f', ('time'))
-                    newvar[:] = tvals
-                    basicvardef(newvar, 'time', 'Time', urefvals)
-                    newvar.setncattr('calendar','gregorian')
 
             else:
@@ -28809 +28846 @@
         newvardims = vardims + dimsslice
         newvarshape = varshape + shapeslice
-        print 'Lluis newvardims:', newvardims, 'newvarshape:', newvarshape
+        print 'Lluis newvardims:', newvardims, 'newvarshape:', newvarshape, 'varv:', varv.shape
+        print 'Lluis shapes sint:', sint.shape
 
 #        newvarslice = []
@@ -28832 +28870 @@
         #print allslices
 
+        vout = None
+        voutt = None
+        varvt = None
         varvt = varv.transpose()
-
+        print 'vs:', varvt[7, 69, :]
         if len(vshape) == 3 and len(sN.shape) == 4:
-            voutt = fsci.module_scientific.multi_spacewightstats_in3drk3_slc3v4(     \
+            voutt = fsci.module_scientific.multi_spaceweightstats_in3drk3_slc3v4(    \
               varin=varvt, ngridsin=sNt, gridsin=sint, percentages=spt,              \
               di1=vshape[2], di2=vshape[1], di3=vshape[0], ds1=ssh[3],               \
@@ -28841 +28882 @@
             vout = voutt.transpose()
         elif len(vshape) == 3 and len(sN.shape) == 3:
-            voutt = fsci.module_scientific.multi_spacewightstats_in3drk3_slc3v3(     \
+            voutt = fsci.module_scientific.multi_spaceweightstats_in3drk3_slc3v3(    \
               varin=varvt, ngridsin=sNt, gridsin=sint, percentages=spt,              \
               di1=vshape[2], di2=vshape[1], di3=vshape[0], ds1=ssh[2],               \
@@ -28855 +28896 @@
             quit(-1)
 
+        print "Lluis Let's check !!", vout.shape
+        print 'i j k it:', i, j, k ,it, 'N:', sN[i,j,k], '<>', slicesvals
+        for iv in range(sN[i,j,k]):
+            print iv, 'c:', sin[:,iv,i,j,k], 'p:', sgpa[iv,i,j,k], 'v:',             \
+              varv[it,sin[0,iv,i,j,k],sin[1,iv,i,j,k]]
+        print 'min max mean mean2 std median count _______'
+        print vout[:,it,i,j,k]
+
         ist = 0
         for statn in statns: