Changeset 2294 in lmdz_wrf for trunk/tools

Timestamp:

Jan 28, 2019, 5:48:59 PM (6 years ago)

Author:

lfita

Message:

Fixing indices errors...

Location:

trunk/tools

Files:

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

trunk/tools/module_scientific.f90

@@ -5304 +5304 @@
   
         CALL grid_within_polygon(Nvertex, vertexgrid, dxB, dyB, dxB*dyB, xCBvals, yCBvals,            &
-          NBvertexmax, xBBvals, yBBvals, Ngridsin(ix,iy), gridsin(ix,iy,:,:))
+          NBvertexmax, xBBvals, yBBvals, Ngridsin(ix,iy), gridsin(ix,iy,1:dxyB,:))
    
         IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin)
@@ -5317 +5317 @@
 
         CALL grid_spacepercen_within_reg(Nvertex, vertexgrid, dxB, dyB, NBvertexmax, xBBvals,        &
-          yBBvals, Ngridsin(ix,iy), poinsin, strict, pareas, percentages(ix,iy,:))
+          yBBvals, Ngridsin(ix,iy), poinsin, strict, pareas, percentages(ix,iy,1:dxyB))
 
         ! Filling areas

trunk/tools/nc_var_tools.py

@@ -28225 +28225 @@
                 percens = percenst.transpose()
 
+                # Remember we are in C-like!
+                gridsin = gridsin - 1
+
                 Ngridsinmax = np.max(Ngridsin)
         else:
@@ -28278 +28281 @@
 
                 ainds = np.array(indices)
+                # Remenber we are in C-like ...
+                ainds = ainds - 1
                 if Ngridsin[0,islc] > 0:
-                    
                     gridsin[0,0:Nt,0,islc] = ainds[0,0:Nt]
                     gridsin[1,0:Nt,0,islc] = ainds[1,0:Nt]
@@ -28296 +28300 @@
                 if not gen.searchInlist(onewnc.dimensions, 'fake'): 
                     newdim = onewnc.createDimension('fake',1)
+                    newdim = onewnc.createDimension('slice_fake',1)
+                    newdim = onewnc.createDimension('fakebnds',4)
+                    newvar = onewnc.createVariable('fake','f',('fake'))
+                    newvar[:] = [0.]
+                    basicvardef(newvar, 'fake', 'fake variable for transforming ' +  \
+                      '1D variable without bounds', '-')
+                    newvar = onewnc.createVariable('slice_fake','f',('slice_fake'))
+                    newvar[:] = [0.]
+                    basicvardef(newvar, 'slice_fake', 'slicing along fake variable '+\
+                      'for transforming variable without bounds', '-')
+                    newvar = onewnc.createVariable('slice_fake_bnds','f',            \
+                      ('fakebnds','fake'))
+                    fakebnds = np.zeros((4,1), dtype=np.float)
+                    fakebnds[0,0] = -1
+                    fakebnds[1,0] = 1
+                    fakebnds[2,0] = 1
+                    fakebnds[3,0] = -1
+                    newvar[:] = fakebnds
+                    basicvardef(newvar, 'slice_fake', 'slicing along fake variable '+\
+                      'for transforming variable without bounds', '-')
+                    onewnc.sync()
 
         if not gen.searchInlist(newslcvarns,dn):
@@ -28345 +28370 @@
             anewvar[:] = areas[:]
 
+            aanewvar = onewnc.createVariable(dn+'area','f',tuple(Srgrid),             \
+              fill_value=gen.fillValueF)
+            basicvardef(aanewvar ,dn+'area', "area of the slice", vu)
+            aanewvar.setncattr('coordinates',' '.join(Srgrid[::-1]))
+
             pnewvar = onewnc.createVariable(dn+'gridpercen','f',tuple(Spgrid),       \
               fill_value=gen.fillValueF)
@@ -28357 +28387 @@
                           gridsin[:,0:Ngridsin[j,i],j,i]
                         pnewvar[0:Ngridsin[j,i],i]= percens[0:Ngridsin[j,i],j,i]
+                        slicearea = 0.
+                        for iv in range(Ngridsin[j,i]):
+                            ix = gridsin[1,iv,i]
+                            iy = gridsin[0,iv,i]
+                            slicearea = slicearea + areas[iy,ix]*percens[iv,j,i]
+                        aanewvar[i]= slicearea
             elif len(dv.shape) == 2:
                 for j in range(dyref):
                    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))
                         innewvar[:,0:Ngridsin[j,i],j,i] =                            \
                           gridsin[:,0:Ngridsin[j,i],j,i]
                         pnewvar[0:Ngridsin[j,i],j,i]= percens[0:Ngridsin[j,i],j,i]
-    
+                        slicearea = 0.
+                        aa = 0.
+                        pp = 0.
+                        for iv in range(Ngridsin[j,i]):
+                            ix = gridsin[1,iv,j,i]
+                            iy = gridsin[0,iv,j,i]
+                            apolygon[:,0] = xdimvarbnds2D[:,iy,ix]
+                            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]
+                            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()
+
             onewnc.sync()
 
@@ -28456 +28519 @@
                     Nin = NpointsAB[jB,iB,jA,iA]
                     innewvar[:,0:Nin,jB,iB,jA,iA] = pointsAB[:,0:Nin,jB,iB,jA,iA]
-                    print Nin, 'Lluis shapes sliceaA:', osliceaA.shape, 'sliceaB:', osliceaB.shape, \
-                    'pnewvar:', pnewvar.shape, 'anewvar:', anewvar.shape
                     for iv in range(Nin):
-                        print 'Lluis jB, iB, jA, iA:', jB,iB,jA,iA
-                        print 'Lluis shapes inpointsA:', inpointsA.shape, 'inpointsB:', inpointsB.shape
-                        print 'Lluis inpointsA:', inpointsA[iv,jA,iA], 'inpointsB:', inpointsB[iv,jB,iB]
                         pA = oslicepA[inpointsA[iv,jA,iA],jA,iA]
                         pB = oslicepB[inpointsB[iv,jB,iB],jB,iB]
-                        print 'Lluis shapes pnewvar:', pnewvar.shape, 'pA:', pA, 'pB:', pB
                         pnewvar[iv,jB,iB,jA,iA]= pA*pB
                         ixA = sliceinA[1,inpointsA[iv,jA,iA],jA,iA]
                         iyA = sliceinA[0,inpointsA[iv,jA,iA],jA,iA]
-                        print '  A:', ixA, iyA
                         aA = osliceaA[iyA,ixA]
                         ixB = sliceinB[1,inpointsB[iv,jB,iB],jB,iB]
                         iyB = sliceinB[0,inpointsB[iv,jB,iB],jB,iB]
-                        print '  B:', ixB, iyB
                         aB = osliceaB[iyB,ixB]
                         anewvar[jB,iB,jA,iA]= (aA*pA)*(aB*pB)