Changeset 2300 in lmdz_wrf

Timestamp:

Jan 29, 2019, 10:50:37 PM (6 years ago)

Author:

lfita

Message:

Trying to underestand why is not working !?

Location:

trunk/tools

Files:

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

trunk/tools/module_scientific.f90 ¶

@@ -6274 +6274 @@
 
     ! Local
-    INTEGER                                              :: ixA, iyA, ixB, iyB, iv
+    INTEGER                                              :: ixA, iyA, ixB, iyB, iv, ii
     INTEGER                                              :: NA, NB
     INTEGER, DIMENSION(dxyA)                             :: ptsA, ptsB
@@ -6313 +6313 @@
               inpointsA(ixA,iyA,iv) = ptsA(iv)
               inpointsB(ixB,iyB,iv) = ptsB(iv)
+              IF (Npoints(ixA,iyA,ixB,iyB) == 4) THEN
+                PRINT *,ixA, iyA,ixB,iyB,'________________'
+                DO ii = 1, Npoints(ixA,iyA,ixB,iyB)
+                  PRINT *,'Lluis ii:', ii, points(ixA,iyA,ixB,iyB,ii,:)
+                  PRINT *,'Lluis A:', inpointsA(ixA,iyA,iv), pointsA(ixA,iyA,ptsA(iv),:)
+                  PRINT *,'Lluis B:', inpointsB(ixB,iyB,iv), pointsB(ixB,iyB,ptsB(iv),:)
+                END DO
+              END IF
             END DO
           END DO

trunk/tools/nc_var_tools.py ¶

@@ -28309 +28309 @@
                 elif (len(unma.shape) == 2):
                     drng = unma.shape[0]*unma.shape[1]
-                    Nt,indices=fsci.module_scientific.multi_index_mat2di(d12=drng,   \
+                    Nt,indicest=fsci.module_scientific.multi_index_mat2di(d12=drng,  \
                       mat=unmat, value=1, d1=unma.shape[1], d2=unma.shape[0])
+                    indices = indicest.transpose()
                 elif (len(unma.shape) == 3):
                     drng = unma.shape[0]*unma.shape[1]*unma.shape[2]
-                    Nt,indices=fsci.module_scientific.multi_index_mat3di(d123=drng,  \
+                    Nt,indicest=fsci.module_scientific.multi_index_mat3di(d123=drng, \
                       mat=unmat, value=1, d1=unma.shape[2], d2=unma.shape[1],        \
                       d3=unma.shape[0])
+                    indices = indicest.transpose()
                 elif (len(unma.shape) == 4):
                     drng = unma.shape[0]*unma.shape[1]*unma.shape[2]*unma.shape[3]
-                    Nt,indices=fsci.module_scientific.multi_index_mat4di(d1234=drng, \
+                    Nt,indicest=fsci.module_scientific.multi_index_mat4di(d1234=drng,\
                       mat=unmat, value=1, d1=unma.shape[3], d2=unma.shape[2],        \
                       d3=unma.shape[1], d4=unma.shape[0])
+                    indices = indicest.transpose()
                 else:
                   rS = str(len(unma.shape))
@@ -28331 +28334 @@
                 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]
+                if Nt > 0:
+                    print 'Lluis:', islc, ':', Nt
+                    gridsin[0,0:Nt,0,islc] = ainds[0:Nt,0]
+                    gridsin[1,0:Nt,0,islc] = ainds[0:Nt,1]
                     percens[0:Nt,0,islc] = 1.
                     for iv in range(Nt):
-                        i = ainds[0,iv]
-                        j = ainds[1,iv]
+                        i = ainds[iv,0]
+                        j = ainds[iv,1]
                         areas[j,i] = 1.
+                        if j == 94: print islc,'Lluis:', j,i,areas[j,i]
 
             # Values for file
@@ -28422 +28427 @@
             anewvar.setncattr('coordinates',' '.join(varfinaldims[::-1]))
             anewvar[:] = areas[:]
+            if dn == 'HGT': 
+                for i in range(areas.shape[1]):
+                    if areas[94,i] == 1.: print 'Lluis 2:', 94,i,areas[94,i]
 
             aanewvar = onewnc.createVariable(dn+'area','f',tuple(Srgrid),             \
@@ -28525 +28533 @@
     sliceinBt = sliceinB.transpose()
 
+    print 'Lluis shapes sliceNAt:', sliceNAt.shape, 'sliceinAt:', sliceinAt.shape
+    print 'Lluis shapes sliceNBt:', sliceNBt.shape, 'sliceinBt:', sliceinBt.shape
+
     NpointsABt, pointsABt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \
       npointsa=sliceNAt, pointsa=sliceinAt, npointsb=sliceNBt, pointsb=sliceinBt,    \
@@ -28534 +28545 @@
 
     # Remembering that it is python (C-like...)
-    inpointsA = inpointsA-1
-    inpointsB = inpointsB-1
+#    inpointsA = inpointsA
+#    inpointsB = inpointsB
+
+    jB = 0
+    iB = 3
+    jA = 10
+    iA = 2
+    print 'iA: _____ '
+    for iv in range(sliceNA[jA,iA]):
+        print iv, sliceinA[:,iv,jA,iA]
+    print 'iB: _____ '
+    for iv in range(sliceNB[jB,iB]):
+        print iv, sliceinB[:,iv,jB,iB]
+    for iv in range(sliceNA[jA,iA]):
+        for iv2 in range(sliceNB[jB,iB]):
+            if sliceinA[0,iv,jA,iA] == sliceinB[0,iv2,jB,iB] and sliceinA[1,iv,jA,iA] == sliceinB[1,iv2,jB,iB]:
+                print iv,':', sliceinA[:,iv,jA,iA], ';', iv2, ':', sliceinB[:,iv2,jB,iB]
+
+    print 'Lluis jB iB jA iA _______',jB,iB,jA,iA,':',NpointsAB[jB,iB,jA,iA] 
+    print 'Lluis inpA shapes', inpointsA.shape, ':', inpointsA[:,jA,iA]
+    print 'Lluis inpB shapes', inpointsB.shape, ':', inpointsB[:,jB,iB]
+    for iv in range(NpointsAB[jB,iB,jA,iA]):
+        print iv, 'NA:', sliceNA[jA,iA], 'pA:', inpointsA[iv,jA,iA], 'c:', sliceinA[:,inpointsA[iv,jA,iA],jA,iA], \
+          'NB:', sliceNB[jB,iB], 'pB:', inpointsB[iv,jB,iB], 'c:', sliceinB[:,inpointsB[iv,jB,iB],jB,iB], '=',    \
+          pointsAB[:,iv,jB,iB,jA,iA]
+
+    quit()
 
     maxNpointsAB = np.max(NpointsAB)
@@ -28579 +28615 @@
                 for iB in range(dxB):
                     Nin = NpointsAB[jB,iB,jA,iA]
+                    print 'jB iB jA iA:', jB,iB,jA,iA, ':', Nin
                     innewvar[:,0:Nin,jB,iB,jA,iA] = pointsAB[:,0:Nin,jB,iB,jA,iA]
                     slicearea = 0.
@@ -28592 +28629 @@
                         aB = osliceaB[iyB,ixB]
                         # I do not understand why this is needed !!
-                        if aA != gen.mamat[1] and aB != gen.mamat[1]: 
-                            anewvar[iv.jB,iB,jA,iA]= (aA*pA)*(aB*pB)
+                        if type(aA)!=type(gen.mamat[1]) and                          \
+                          type(aB)!=type(gen.mamat[1]):
+                            anewvar[iv,jB,iB,jA,iA]= (aA*pA)*(aB*pB)
                             slicearea = slicearea + (aA*pA)*(aB*pB)
+                        else:
+                            print type(aA) != type(gen.mamat[1]), type(aB) != type(gen.mamat[1])
+                            print 'Lluis aA ________', iyA, ixA, 'shape osliceaA:', osliceaA.shape, ':', aA, gen.mamat[1]
+                            print osliceaA[iyA-2:iyA+3,ixA-2:ixA+3]
+                            print 'Lluis aB ________', iyB, ixB, 'shape osliceaB:', osliceaB.shape, ':', aB
+                            print osliceaB[iyB-2:iyB+3,ixB-2:ixB+3]
+                            #quit()
                         aanewvar[jB,iB,jA,iA]= slicearea
 
@@ -28701 +28746 @@
                             aB = osliceaB[iyB,ixB]
                             # I do not understand why this is needed !!
-                            if aA != gen.mamat[1] and aB != gen.mamat[1]: 
-                                anewvar[iv,jB,iB,jA,iA]= (aA*pA)*(aB*pB)
+                            if type(aA)!=type(gen.mamat[1]) and                          \
+                              type(aB)!=type(gen.mamat[1]):
                                 slicearea = slicearea + (aA*pA)*(aB*pB)
                             aanewvar[jB,iB,jA,iA]= slicearea
@@ -28731 +28776 @@
     # Percentages as the areal fraction of each grid with respect slice size
     sgp = sga/sa
+    i=3
+    j=10
+    k=2
+    print 'Lluis percentages ________', i, j,k
+    Npt = sN[i,j,k]
+    for iv in range(Npt):
+        print 'Lluis:', iv,'pt:', sin[:,iv,i,j,k], 'a:', sa[i,j,k], 'ga:', sga[iv,i,j,k]
+    quit()
 
     ssh = sN.shape
@@ -28841 +28894 @@
             newvar = onewvars[statn]
             print 'Lluis stn:', statn, 'shape:', newvar.shape
-            newvar[:] = vout[ist]
+            rightvals = np.where(vout[6] == 0, gen.fillValueF, vout[ist])
+            newvar[:] = rightvals
             onewnc.sync()
+            ist = ist + 1
 
         quit()