Changeset 2311 in lmdz_wrf


Ignore:
Timestamp:
Feb 1, 2019, 7:59:22 PM (6 years ago)
Author:
lfita
Message:

Fixing:

  • NaN test on area of the slices
  • memory overload
  • 1-value statistics weighted computations
Location:
trunk/tools
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/tools/module_scientific.f90

    r2308 r2311  
    57015701        DO s3 =1, ds3
    57025702          Nin = Ngridsin(s1,s2,s3)
    5703           IF (ALLOCATED(gin)) DEALLOCATE(gin)
    5704           ALLOCATE(gin(Nin,2))
    5705           IF (ALLOCATED(pin)) DEALLOCATE(pin)
    5706           ALLOCATE(pin(Nin))
    5707           IF (ALLOCATED(vin)) DEALLOCATE(vin)
    5708           ALLOCATE(vin(Nin,di3))
    5709           IF (ALLOCATED(svin)) DEALLOCATE(svin)
    5710           ALLOCATE(svin(Nin))
    5711           gin = gridsin(s1,s2,s3,1:Nin,:)
    5712           pin = percentages(s1,s2,s3,1:Nin)
    5713 
    5714           ! Getting the values
    5715           DO iv=1, Nin
    5716             i1 = gin(iv,1)
    5717             i2 = gin(iv,2)
    5718             vin(iv,:) = varin(i1,i2,:)
    5719           END DO
    5720 
    57215703          ! Computing along d3
    5722           DO i3=1, di3
    5723             minv = fillVal64
    5724             maxv = -fillVal64
    5725             meanv = zeroRK
    5726             mean2v = zeroRK
    5727             stdv = zeroRK
    5728             minv = MINVAL(vin(:,i3))
    5729             maxv = MAXVAL(vin(:,i3))
    5730             meanv = SUM(vin(:,i3)*pin)
    5731             mean2v = SUM(vin(:,i3)**2*pin)
    5732             DO iv=1,Nin
    5733               stdv = stdv + ( (meanv - vin(iv,i3))*pin(iv) )**2
     5704          IF (Nin > 1) THEN
     5705            IF (ALLOCATED(gin)) DEALLOCATE(gin)
     5706            ALLOCATE(gin(Nin,2))
     5707            IF (ALLOCATED(pin)) DEALLOCATE(pin)
     5708            ALLOCATE(pin(Nin))
     5709            IF (ALLOCATED(vin)) DEALLOCATE(vin)
     5710            ALLOCATE(vin(Nin,di3))
     5711            IF (ALLOCATED(svin)) DEALLOCATE(svin)
     5712            ALLOCATE(svin(Nin))
     5713            gin = gridsin(s1,s2,s3,1:Nin,:)
     5714            pin = percentages(s1,s2,s3,1:Nin)
     5715
     5716            ! Getting the values
     5717            DO iv=1, Nin
     5718              i1 = gin(iv,1)
     5719              i2 = gin(iv,2)
     5720              vin(iv,:) = varin(i1,i2,:)
    57345721            END DO
    5735             stdv = SQRT(stdv)
    5736             svin = vin(:,i3)
    5737             CALL SortR_K(svin, Nin)
    5738             medv = svin(INT(Nin/2))
    5739             varout(s1,s2,s3,i3,1) = minv
    5740             varout(s1,s2,s3,i3,2) = maxv
    5741             varout(s1,s2,s3,i3,3) = meanv
    5742             varout(s1,s2,s3,i3,4) = mean2v
    5743             varout(s1,s2,s3,i3,5) = stdv
    5744             varout(s1,s2,s3,i3,6) = medv
    5745             varout(s1,s2,s3,i3,7) = Nin*1.
    5746           END DO
     5722            DO i3=1, di3
     5723              minv = fillVal64
     5724              maxv = -fillVal64
     5725              meanv = zeroRK
     5726              mean2v = zeroRK
     5727              stdv = zeroRK
     5728              minv = MINVAL(vin(:,i3))
     5729              maxv = MAXVAL(vin(:,i3))
     5730              meanv = SUM(vin(:,i3)*pin)
     5731              mean2v = SUM(vin(:,i3)**2*pin)
     5732              DO iv=1,Nin
     5733                stdv = stdv + ( (meanv - vin(iv,i3))*pin(iv) )**2
     5734              END DO
     5735              stdv = SQRT(stdv)
     5736              svin = vin(:,i3)
     5737              CALL SortR_K(svin, Nin)
     5738              medv = svin(INT(Nin/2))
     5739              varout(s1,s2,s3,i3,1) = minv
     5740              varout(s1,s2,s3,i3,2) = maxv
     5741              varout(s1,s2,s3,i3,3) = meanv
     5742              varout(s1,s2,s3,i3,4) = mean2v
     5743              varout(s1,s2,s3,i3,5) = stdv
     5744              varout(s1,s2,s3,i3,6) = medv
     5745              varout(s1,s2,s3,i3,7) = Nin*1.
     5746            END DO
     5747          ELSE
     5748            varout(s1,s2,s3,:,1) = varin(i1,i2,:)
     5749            varout(s1,s2,s3,:,2) = varin(i1,i2,:)
     5750            varout(s1,s2,s3,:,3) = varin(i1,i2,:)
     5751            varout(s1,s2,s3,:,4) = varin(i1,i2,:)*varin(i1,i2,:)
     5752            varout(s1,s2,s3,:,5) = zeroRK
     5753            varout(s1,s2,s3,:,6) = varin(i1,i2,:)
     5754            varout(s1,s2,s3,:,7) = Nin*1.
     5755          END IF
    57475756        END DO
    57485757      END DO
     
    58175826          DO s4 =1, ds4
    58185827            Nin = Ngridsin(s1,s2,s3,s4)
    5819             IF (ALLOCATED(gin)) DEALLOCATE(gin)
    5820             ALLOCATE(gin(Nin,2))
    5821             IF (ALLOCATED(pin)) DEALLOCATE(pin)
    5822             ALLOCATE(pin(Nin))
    5823             IF (ALLOCATED(vin)) DEALLOCATE(vin)
    5824             ALLOCATE(vin(Nin,di3))
    5825             IF (ALLOCATED(svin)) DEALLOCATE(svin)
    5826             ALLOCATE(svin(Nin))
    5827             gin = gridsin(s1,s2,s3,s4,1:Nin,:)
    5828             pin = percentages(s1,s2,s3,s4,1:Nin)
    5829 
    5830             minv = fillVal64
    5831             maxv = -fillVal64
    5832             meanv = zeroRK
    5833             mean2v = zeroRK
    5834             stdv = zeroRK
    5835 
    5836             ! Getting the values
    5837             DO iv=1, Nin
    5838               i1 = gin(iv,1)
    5839               i2 = gin(iv,2)
    5840               vin(iv,:) = varin(i1,i2,:)
    5841             END DO
    5842 
    5843             ! Computing along d3
    5844             DO i3=1, di3
    5845               minv = fillVal64
    5846               maxv = -fillVal64
    5847               meanv = zeroRK
    5848               mean2v = zeroRK
    5849               stdv = zeroRK
    5850 
    5851               minv = MINVAL(vin(:,i3))
    5852               maxv = MAXVAL(vin(:,i3))
    5853               meanv = SUM(vin(:,i3)*pin)
    5854               mean2v = SUM(vin(:,i3)**2*pin)
    5855               DO iv=1,Nin
    5856                 stdv = stdv + ( (meanv - vin(iv,i3))*pin(iv) )**2
     5828            IF (Nin > 1) THEN
     5829              IF (ALLOCATED(gin)) DEALLOCATE(gin)
     5830              ALLOCATE(gin(Nin,2))
     5831              IF (ALLOCATED(pin)) DEALLOCATE(pin)
     5832              ALLOCATE(pin(Nin))
     5833              IF (ALLOCATED(vin)) DEALLOCATE(vin)
     5834              ALLOCATE(vin(Nin,di3))
     5835              IF (ALLOCATED(svin)) DEALLOCATE(svin)
     5836              ALLOCATE(svin(Nin))
     5837              gin = gridsin(s1,s2,s3,s4,1:Nin,:)
     5838              pin = percentages(s1,s2,s3,s4,1:Nin)
     5839
     5840              ! Getting the values
     5841              DO iv=1, Nin
     5842                i1 = gin(iv,1)
     5843                i2 = gin(iv,2)
     5844                vin(iv,:) = varin(i1,i2,:)
    58575845              END DO
    5858               stdv = SQRT(stdv)
    5859               svin = vin(:,i3)
    5860               CALL SortR_K(svin, Nin)
    5861               medv = svin(INT(Nin/2))
    5862               varout(s1,s2,s3,s4,i3,1) = minv
    5863               varout(s1,s2,s3,s4,i3,2) = maxv
    5864               varout(s1,s2,s3,s4,i3,3) = meanv
    5865               varout(s1,s2,s3,s4,i3,4) = mean2v
    5866               varout(s1,s2,s3,s4,i3,5) = stdv
    5867               varout(s1,s2,s3,s4,i3,6) = medv
    5868               varout(s1,s2,s3,s4,i3,7) = Nin*1.
    5869             END DO
     5846              ! Computing along d3
     5847              DO i3=1, di3
     5848                minv = fillVal64
     5849                maxv = -fillVal64
     5850                meanv = zeroRK
     5851                mean2v = zeroRK
     5852                stdv = zeroRK
     5853
     5854                minv = MINVAL(vin(:,i3))
     5855                maxv = MAXVAL(vin(:,i3))
     5856                meanv = SUM(vin(:,i3)*pin)
     5857                mean2v = SUM(vin(:,i3)**2*pin) 
     5858                DO iv=1,Nin
     5859                  stdv = stdv + ( (meanv - vin(iv,i3))*pin(iv) )**2
     5860                END DO
     5861                stdv = SQRT(stdv)
     5862                svin = vin(:,i3)
     5863                CALL SortR_K(svin, Nin)
     5864                medv = svin(INT(Nin/2))
     5865                varout(s1,s2,s3,s4,i3,1) = minv
     5866                varout(s1,s2,s3,s4,i3,2) = maxv
     5867                varout(s1,s2,s3,s4,i3,3) = meanv
     5868                varout(s1,s2,s3,s4,i3,4) = mean2v
     5869                varout(s1,s2,s3,s4,i3,5) = stdv
     5870                varout(s1,s2,s3,s4,i3,6) = medv
     5871                varout(s1,s2,s3,s4,i3,7) = Nin*1.
     5872              END DO
     5873            ELSE
     5874                varout(s1,s2,s3,s4,:,1) = varin(i1,i2,:)
     5875                varout(s1,s2,s3,s4,:,2) = varin(i1,i2,:)
     5876                varout(s1,s2,s3,s4,:,3) = varin(i1,i2,:)
     5877                varout(s1,s2,s3,s4,:,4) = varin(i1,i2,:)*varin(i1,i2,:)
     5878                varout(s1,s2,s3,s4,:,5) = zeroRK
     5879                varout(s1,s2,s3,s4,:,6) = varin(i1,i2,:)
     5880                varout(s1,s2,s3,s4,:,7) = Nin*1.
     5881            END IF
    58705882          END DO
    58715883        END DO
  • trunk/tools/nc_var_tools.py

    r2308 r2311  
    2848828488                            aa = aa + areas[iy,ix]
    2848928489                            pp = pp + percens[iv,j,i]
     28490                        if np.isnan(slicearea):
     28491                            print errormsg
     28492                            print '  ' + fname + ': Nan slice area for slice:',          \
     28493                              j, i, 'N grids:', Ngridsin[j,i],' !!'
     28494                            print '     #grid iy,ix grid_area grid_percen______'
     28495                            for iv in range(Ngridsin[j,i]):
     28496                                ix = gridsin[1,iv,j,i]
     28497                                iy = gridsin[0,iv,j,i]
     28498                                print iv, iy, ix, areas[iy,ix], percens[j,i]
     28499                            quit(-1)
    2849028500                        aanewvar[j,i]= slicearea
    2849128501                        # This should be true for the regular lon/lat projections
     
    2862328633                        anewvar[iv,jB,iB,jA,iA]= (aA*pA)*(aB*pB)
    2862428634                        slicearea = slicearea + (aA*pA)*(aB*pB)
     28635                    if np.isnan(slicearea):
     28636                        print errormsg
     28637                        print '  ' + fname + ': Nan slice area for slice:',          \
     28638                          jA, iA, jB, iB, 'N grids:', Nin,' !!'
     28639                        print '     #grid A grid_area A grid_percen B grid_area ' +  \
     28640                          'B grid_percen______'
     28641                        for iv in range(Nin):
     28642                            print iv, aA, pA, aB, pB
     28643                        quit(-1)
    2862528644                    aanewvar[jB,iB,jA,iA]= slicearea
    2862628645
     
    2878628805                              ' slice:', jA, iA, jB, iB, 'N grids:', Nin,' !!'
    2878728806                            print '     slice total area:', aanewvar[jA,iA,jB,iB]
    28788                             print '     #grid grid_area grid_percen _______'
     28807                            print '     #grid grid_area slice_area grid_percen ______'
    2878928808                            for iv in range(Nin):
    28790                                 print iv, anewvar[iv,jA,iA,jB,iB], sgpa[iv,jA,iA,jB,iB]
     28809                                print iv, anewvar[iv,jA,iA,jB,iB],                   \
     28810                                  aanewvar[jA,iA,jB,iB], sgpa[iv,jA,iA,jB,iB]
    2879128811                            quit(-1)
    2879228812    onewnc.sync()
     
    2883328853        print '    ' + varn + ' ...'
    2883428854        ovar = onc.variables[varn]
    28835         varv = ovar[:]
    2883628855        vdimns = list(ovar.dimensions)
    2883728856        vshape = list(ovar.shape)
     
    2890328922        voutt = None
    2890428923        varvt = None
    28905         varvt = varv.transpose()
    28906         if len(vshape) == 3 and len(sN.shape) == 4:
    28907             voutt = fsci.module_scientific.multi_spaceweightstats_in3drk3_slc3v4(    \
    28908               varin=varvt, ngridsin=sNt, gridsin=sint, percentages=spt,              \
    28909               di1=vshape[2], di2=vshape[1], di3=vshape[0], ds1=ssh[3],               \
    28910               ds2=ssh[2], ds3=ssh[1], ds4=ssh[0], maxngridsin=sin.shape[1])
    28911             vout = voutt.transpose()
    28912         elif len(vshape) == 3 and len(sN.shape) == 3:
    28913             voutt = fsci.module_scientific.multi_spaceweightstats_in3drk3_slc3v3(    \
    28914               varin=varvt, ngridsin=sNt, gridsin=sint, percentages=spt,              \
    28915               di1=vshape[2], di2=vshape[1], di3=vshape[0], ds1=ssh[2],               \
    28916               ds2=ssh[1], ds3=ssh[0], maxngridsin=sin.shape[1])
    28917             vout = voutt.transpose()
    28918         else:
    28919             print errormsg
    28920             print '  ' + fname + ": variable rank:", len(vshape), 'combined with',   \
    28921              ' slice rank:', len(sN.shape), 'not ready!!'
    28922             print '    available ones _______'
    28923             print ' var_rank: 3   slice_rank: 3'
    28924             print ' var_rank: 3   slice_rank: 4'
    28925             quit(-1)
    28926 
    28927         #print "Lluis Let's check !!", vout.shape
    28928         #print 'i j k it:', i, j, k ,it, 'N:', sN[i,j,k], '<>', slicesvals
    28929         #for iv in range(sN[i,j,k]):
    28930         #    print iv, 'c:', sin[:,iv,i,j,k], 'p:', sgpa[iv,i,j,k], 'v:',             \
    28931         #      varv[it,sin[0,iv,i,j,k],sin[1,iv,i,j,k]]
    28932         #print 'min max mean mean2 std median count _______'
    28933         #print vout[:,it,i,j,k]
    28934 
    28935         ist = 0
    28936         for statn in statns:
    28937             newvar = onewvars[statn]
    28938             rightvals = np.where(vout[6] == 0, gen.fillValueF, vout[ist])
    28939             newvar[:] = rightvals
    28940             onewnc.sync()
    28941             ist = ist + 1
    28942 
    28943         onewnc.sync()
     28924
     28925        # Checking for size...
     28926        Nvdims = len(vshape)
     28927        varslc = {}
     28928        if Nvdims == 3:
     28929            vSIZE = np.prod(varv.shape[:])*8.
     28930            if vSIZE > gen.oneGB*3.:
     28931                print '  ' + infmsg
     28932                print '    variable size: ', vSIZE / gen.oneGB,' larger than 3GB !!'
     28933                print '    splitting time evolution to compute'
     28934                vHSIZE = np.prod(varv.shape[1:3])*8.
     28935                print '    var. horizontal size:', vHSIZE, 'dimt:', varv.shape[0]
     28936                tsteps = gen.oneGB*3. / vHSIZE
     28937                tslices = range(0,varv.shape[2],tsteps)
     28938                print '    time-slices:', tklices
     28939                for itt in tslcices[1,len(tslices)]:
     28940                    iit = tslices[itt-1]
     28941                    eet = tslices[itt]
     28942                    varslc[itt-1] = [slice(iit,eet), slice(0,vshape[1]),             \
     28943                      slice(0,vbshape[2])]
     28944            else:
     28945                tslices = range(0,varv.shape[0])
     28946                varslc[0] = [slice(0,vshape[0]), slice(0,vshape[1]),                 \
     28947                  slice(0,vshape[2])]
     28948
     28949        for itt in varslc.keys():
     28950           slcv = varslc[itt]
     28951           varv = ovar[tuple(slcv)]
     28952           d0 = varv.shape[0]
     28953           d1 = varv.shape[1]
     28954           d2 = varv.shape[2]
     28955           varvt = varv.transpose()
     28956           if Nvdims == 3 and len(sN.shape) == 4:
     28957               voutt=fsci.module_scientific.multi_spaceweightstats_in3drk3_slc3v4(   \
     28958                 varin=varvt, ngridsin=sNt, gridsin=sint, percentages=spt,           \
     28959                 di1=d2, di2=d1, di3=d0, ds1=ssh[3], ds2=ssh[2], ds3=ssh[1],         \
     28960                 ds4=ssh[0], maxngridsin=sin.shape[1])
     28961           elif Nvdims == 3 and len(sN.shape) == 3:
     28962               voutt=fsci.module_scientific.multi_spaceweightstats_in3drk3_slc3v3(   \
     28963                 varin=varvt, ngridsin=sNt, gridsin=sint, percentages=spt,           \
     28964                 di1=d2, di2=d1, di3=d0, ds1=ssh[2], ds2=ssh[1], ds3=ssh[0],         \
     28965                 maxngridsin=sin.shape[1])
     28966           else:
     28967               print errormsg
     28968               print '  ' + fname + ": variable rank:", Nvdims, 'combined with',     \
     28969                 ' slice rank:', len(sN.shape), 'not ready!!'
     28970               print '    available ones _______'
     28971               print ' var_rank: 3   slice_rank: 3'
     28972               print ' var_rank: 3   slice_rank: 4'
     28973               quit(-1)
     28974
     28975           vout = voutt.transpose()
     28976           ist = 0
     28977           for statn in statns:
     28978               newvar = onewvars[statn]
     28979               rightvals = np.where(vout[6] == 0, gen.fillValueF, vout[ist])
     28980               newvar[tuple(slcv)] = rightvals
     28981               onewnc.sync()
     28982               ist = ist + 1
     28983           onewnc.sync()
     28984           #print "Lluis Let's check !!", vout.shape
     28985           #print 'i j k it:', i, j, k ,it, 'N:', sN[i,j,k], '<>', slicesvals
     28986           #for iv in range(sN[i,j,k]):
     28987           #    print iv, 'c:', sin[:,iv,i,j,k], 'p:', sgpa[iv,i,j,k], 'v:',          \
     28988           #      varv[it,sin[0,iv,i,j,k],sin[1,iv,i,j,k]]
     28989           #print 'min max mean mean2 std median count _______'
     28990           #print vout[:,it,i,j,k]
    2894428991
    2894528992    # Add global attributes
Note: See TracChangeset for help on using the changeset viewer.