Changeset 2295 in lmdz_wrf for trunk/tools

Timestamp:

Jan 28, 2019, 10:49:39 PM (6 years ago)

Author:

lfita

Message:

Incroporation of the slice size when NO bounds

Location:

trunk/tools

Files:

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

trunk/tools/module_scientific.f90

@@ -5263 +5263 @@
 ! Local
    INTEGER                                               :: iv, ix, iy
-   INTEGER                                               :: Nvertex
+   INTEGER                                               :: Nvertex, Nptin
    INTEGER, ALLOCATABLE, DIMENSION(:,:)                  :: poinsin
    CHARACTER(len=20)                                     :: IS
@@ -5316 +5316 @@
         END DO
 
+        Nptin = Ngridsin(ix,iy)
         CALL grid_spacepercen_within_reg(Nvertex, vertexgrid, dxB, dyB, NBvertexmax, xBBvals,        &
-          yBBvals, Ngridsin(ix,iy), poinsin, strict, pareas, percentages(ix,iy,1:dxyB))
+          yBBvals, Ngridsin(ix,iy), poinsin, strict, pareas, percentages(ix,iy,1:Nptin))
 
         ! Filling areas

trunk/tools/nc_var_tools.py

@@ -28231 +28231 @@
         else:
             dn = varn + ''
+            varslcv = slicesinf[dn]
+            Nslices = varslcv[0]
+            ovardims = varslcv[1]
+            slcvar = varslcv[2]
+            lcvalsc = varslcv[3]
+            slcvals = varslcv[4]
+
             dv = np.zeros((1,Nslices), dtype=int)
             # Here is only done by 1D variables !!!
@@ -28236 +28243 @@
             print '  ' + fname + ": slicing variable '" + dn + "' without bounds !!"
             print '    directly using grid point values'
-
-            #if len(varv.shape) > 1: 
-            #    print errormsg
-            #    print '  ' + fname + ": wrong shape for slicing vriable '" + varn + "' !!"
-            #    print '    function only ready to slice 1D varables and it has:',        \
-            #      varv.shape
-            #    quit(-1)
+            print '    remaining dimensions:', rmdims
+
+            refblon1D = np.zeros((Nslices,4), dtype=np.float)
+            for iv in range(Nslices):
+                slcvv = slcvals[iv]
+                refblon1D[iv,0] = slcvv[0]
+                refblon1D[iv,1] = slcvv[1]
+                refblon1D[iv,2] = slcvv[1]
+                refblon1D[iv,3] = slcvv[0]
+            refblat1D = np.zeros((1,4), dtype=np.float)
+            refblat1D[0,:] = np.array([-1,1,1,-1])
+
+            # Getting boundaries to compute area of the slices
+            iid = 0
+            for dimn in rmdims:
+                if slicebndsdim.has_key(dimn): 
+                    print "    bounded dimension '" + dimn + "' by '" + slicebndsdim[dimn] + "' !!"
+                    if iid == 0:
+                        ogetblat = onc.variables[slicebndsdim[dimn]]
+                        # Removing undesired dimensions from bounds slicing variable
+                        if sliceremovedim is not None:
+                            getblat, rmd, rms = ovar_reducedims(ogetblat, sliceremovedim)
+                        else:
+                            getblat = ogetblat[:]
+                            rmd = ogetblat.dimensions
+                            rms = list(ogetblat.shape)
+                    elif iid == 1:
+                        ogetblon = onc.variables[slicebndsdim[dimn]]
+                        # Removing undesired dimensions from bounds slicing variable
+                        if sliceremovedim is not None:
+                            getblon, rmd, rms = ovar_reducedims(ogetblon, sliceremovedim)
+                        else:
+                            getblon = ogetblon[:]
+                            rmd = ogetblon.dimensions
+                            rms = list(ogetblon.shape)
+                    iid = iid + 1
+            xdimvarbnds2D = getblon
+            ydimvarbnds2D = getblat
 
             Ngridsin = np.zeros((1,Nslices), dtype=int)           
@@ -28253 +28291 @@
             areas = np.ones(tuple(varfinalshape), dtype=np.float)
             percens = np.zeros((Ngridsinmax,1,Nslices), dtype=np.float)
+            
             for islc in range(Nslices):
                 unmasked = slcvar[islc,]
@@ -28283 +28322 @@
                 # Remenber we are in C-like ...
                 ainds = ainds - 1
+
                 if Ngridsin[0,islc] > 0:
                     gridsin[0,0:Nt,0,islc] = ainds[0,0:Nt]
@@ -28397 +28437 @@
                    for i in range(dxref):
                         apolygon = np.zeros((4,2), dtype=np.float)
-                        apolygon[0,0] = refblon1D[i,0]
-                        apolygon[1,0] = refblon1D[i,1]
-                        apolygon[2,0] = refblon1D[i,1]
-                        apolygon[3,0] = refblon1D[i,0]
-                        apolygon[0,1] = refblat1D[j,1]
-                        apolygon[1,1] = refblat1D[j,1]
-                        apolygon[2,1] = refblat1D[j,0]
-                        apolygon[3,1] = refblat1D[j,0]
-                        ap = np.abs(fsci.module_scientific.shoelace_area_polygon(nvertex=4, poly=apolygon))
+                        #apolygon[0,0] = refblon1D[i,0]
+                        #apolygon[1,0] = refblon1D[i,1]
+                        #apolygon[2,0] = refblon1D[i,1]
+                        #apolygon[3,0] = refblon1D[i,0]
+                        #apolygon[0,1] = refblat1D[j,1]
+                        #apolygon[1,1] = refblat1D[j,1]
+                        #apolygon[2,1] = refblat1D[j,0]
+                        #apolygon[3,1] = refblat1D[j,0]
+                        #ap = np.abs(fsci.module_scientific.shoelace_area_polygon(nvertex=4, poly=apolygon))
                         innewvar[:,0:Ngridsin[j,i],j,i] =                            \
                           gridsin[:,0:Ngridsin[j,i],j,i]
@@ -28418 +28458 @@
                             apolygon[:,1] = ydimvarbnds2D[:,iy,ix]
                             ag = np.abs(fsci.module_scientific.shoelace_area_polygon(nvertex=4, poly=apolygon))
-                            print '  lluis iv:', iv, ':', iy, ix, 'ag:', ag, 'areas:', areas[iy,ix], 'percens:', percens[iv,j,i]
-                            slicearea = slicearea + areas[iy,ix]*percens[iv,j,i]
+                            if varn == 'HGT':
+                                print 'Lluis ref dy:', dyref, 'dx:', dxref 
+                                print '  lluis iv:', iv, ':', iy, ix, '<>', apolygon[1:3,0],'x',apolygon[2:4,1],'ag:', ag, 'areas:', areas[iy,ix], 'percens:', percens[iv,j,i]
+                            slicearea = slicearea + ag*percens[iv,j,i]
                             aa = aa + areas[iy,ix]
                             pp = pp + percens[iv,j,i]
                         aanewvar[j,i]= slicearea
-                        print '  areas Lluis =', refblon1D[i,:], 'x', refblat1D[j,:], \
-                          'j i:', j,i, 'ap', ap, 'aanewvar:', slicearea, 'aa:', aa, 'pp:', pp
-                        if slicearea > ap: quit()
+                        #print '  areas Lluis =', refblon1D[i,:], 'x', refblat1D[j,:], \
+                        #  'j i:', j,i, 'ap', ap, 'aanewvar:', slicearea, 'aa:', aa, 'pp:', pp
+                        if varn == 'HGT':
+                            print '  areas Lluis =', refblon1D[i,:], 'x', refblat1D[j,:], \
+                              'j i:', j,i, 'aanewvar:', slicearea, 'aa:', aa, 'pp:', pp
+                        #if slicearea > ap: quit()
 
             onewnc.sync()