Changeset 1815 in lmdz_wrf for trunk/tools/nc_var_tools.py


Ignore:
Timestamp:
Mar 19, 2018, 6:36:13 PM (7 years ago)
Author:
lfita
Message:

introductino of pinterp on y-axis pieces

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/nc_var_tools.py

    r1811 r1815  
    180180else:
    181181    isFR64 = False
    182 print 'Lluis oneRK:', type(fdef.module_definitions.onerk), 'isFR64:', isFR64
    183182
    184183## Constants
    … …  
    1808518084    gen.check_arguments(fname, values, arguments, ',')
    1808618085
    18087     interplevs = np.array(values.split(',')[0].split(':'), dtype=np.float64)
     18086    if isFR64:
     18087        interplevs = np.array(values.split(',')[0].split(':'), dtype=np.float64)
     18088    else:
     18089        interplevs = np.array(values.split(',')[0].split(':'), dtype=np.float)
    1808818090    linloginterp = np.int32(values.split(',')[1])
    1808918091    extrap = np.int32(values.split(',')[2])
    … …  
    1819118193            pres = pres0.astype('float64')
    1819218194        else:
    18193             pres = pres0.astype('float')
     18195            pres = pres0.copy()
    1819418196
    1819518197        dimx = pres.shape[3]
    … …  
    1820618208            pres = ovar1[:].astype('float64')
    1820718209        else:
    18208             pres = ovar1[:].astype('float')
     18210            pres = ovar1[:]
    1820918211        dimx = pres.shape[3]
    1821018212        dimy = pres.shape[2]
    … …  
    1821918221            pres = ovar1[:].astype('float64')
    1822018222        else:
    18221             pres = ovar1[:].astype('float')
     18223            pres = ovar1[:]
    1822218224        dimx = pres.shape[3]
    1822318225        dimy = pres.shape[2]
    … …  
    1823218234            pres = ovar1[:].astype('float64')
    1823318235        else:
    18234             pres = ovar1[:].astype('float')
     18236            pres = ovar1[:]
    1823518237        dimx = pres.shape[3]
    1823618238        dimy = pres.shape[2]
    … …  
    1824918251            psfc = ovar1[:].astype('float64')
    1825018252        else:
    18251             psfc = ovar1[:].astype('float')
     18253            psfc = ovar1[:]
    1825218254    elif modname == 'LMDZ' or modname == 'cfLMDZ':
    1825318255        ovar1 = onc.variables['psol']
    … …  
    1825518257            psfc = ovar1[:].astype('float64')
    1825618258        else:
    18257             psfc = ovar1[:].astype('float')
     18259            psfc = ovar1[:]
    1825818260    else:
    1825918261        ovar1 = onc.variables['ps']
    … …  
    1826118263            psfc = ovar1[:].astype('float64')
    1826218264        else:
    18263             psfc = ovar1[:].astype('float')
     18265            psfc = ovar1[:]
    1826418266
    1826518267    # geopotential height
    … …  
    1827018272        unstg = list(geop0.shape)
    1827118273        unstg[1] = unstg[1] - 1
    18272         geop = np.zeros(tuple(unstg), dtype=np.float)
    1827318274        if isFR64:
     18275            geop = np.zeros(tuple(unstg), dtype=np.float64)
    1827418276            geop = 0.5*(geop0[:,1:dimz+1,:,:] + geop0[:,0:dimz,:,:]).astype('float64')
    1827518277        else:
    18276             geop = 0.5*(geop0[:,1:dimz+1,:,:] + geop0[:,0:dimz,:,:]).astype('float')
     18278            geop = np.zeros(tuple(unstg), dtype=np.float)
     18279            geop = 0.5*(geop0[:,1:dimz+1,:,:] + geop0[:,0:dimz,:,:])
    1827718280    elif modname == 'LMDZ' or modname == 'cfLMDZ':
    1827818281        ovar1 = onc.variables['geop']
    … …  
    1828018283            geop = ovar1[:].astype('float64')
    1828118284        else:
    18282             geop = ovar1[:].astype('float')
     18285            geop = ovar1[:]
    1828318286    else:
    1828418287        ovar1 = onc.variables['z']
    … …  
    1828618289            geop = ovar1[:].astype('float64')
    1828718290        else:
    18288             geop = ovar1[:].astype('float')
     18291            geop = ovar1[:]
    1828918292
    1829018293    # terrain height
    … …  
    1829418297            hgt = ovar1[0,:,:].astype('float64')
    1829518298        else:
    18296             hgt = ovar1[0,:,:].astype('float')
     18299            hgt = ovar1[0,:,:]
    1829718300    elif modname == 'LMDZ' or modname == 'cfLMDZ':
    1829818301        grav = 9.81
    … …  
    1830118304            hgt = (ovar1[0,:,:]/grav).astype('float64')
    1830218305        else:
    18303             hgt = (ovar1[0,:,:]/grav).astype('float')
     18306            hgt = (ovar1[0,:,:]/grav)
    1830418307    else:
    1830518308        ovar1 = onc.variables['orog']
    … …  
    1830718310            hgt = ovar1[:].astype('float64')
    1830818311        else:
    18309             hgt = ovar1[:].astype('float')
     18312            hgt = ovar1[:]
    1831018313
    1831118314    # water vapour mixing ratio
    … …  
    1831518318            qv = ovar1[:].astype('float64')
    1831618319        else:
    18317             qv = ovar1[:].astype('float')
     18320            qv = ovar1[:]
    1831818321    elif modname == 'LMDZ' or modname == 'cfLMDZ':
    1831918322        ovar1 = onc.variables['ovap']
    … …  
    1832118324            qv = ovar1[:].astype('float64')
    1832218325        else:
    18323             qv = ovar1[:].astype('float')
     18326            qv = ovar1[:]
    1832418327    else:
    1832518328        ovar1 = onc.variables['hus']
    … …  
    1832718330            qv = ovar1[:].astype('float64')
    1832818331        else:
    18329             qv = ovar1[:].astype('float')
     18332            qv = ovar1[:]
    1833018333
    1833118334    # temperature
    1833218335    if modname == 'WRF':
    18333         Rd = np.float64(287.04)
    18334         Cp = np.float64(7.*Rd/2.)
    18335         RCP = np.float64(Rd/Cp)
    18336         p0 = np.float64(100000.)
     18336        if isFR64:
     18337            Rd = np.float64(287.04)
     18338            Cp = np.float64(7.*Rd/2.)
     18339            RCP = np.float64(Rd/Cp)
     18340            p0 = np.float64(100000.)
     18341        else:
     18342            Rd = 287.04
     18343            Cp = 7.*Rd/2.
     18344            RCP = Rd/Cp
     18345            p0 = 100000.
    1833718346        ovar10 = onc.variables['T']
    1833818347        var10 = ovar10[:]
    … …  
    1834018349            ovar1 = (var10).astype('float64')
    1834118350        else:
    18342             ovar1 = (var10).astype('float')
     18351            ovar1 = var10[:]
    1834318352        temp0 = (ovar1[:]+300.)*(pres[:]/p0)**RCP
    1834418353        if isFR64:
    1834518354            temp = temp0.astype('float64')
    1834618355        else:
    18347             temp = temp0.astype('float')
     18356            temp = temp0.copy()
    1834818357    elif modname == 'LMDZ' or modname == 'cfLMDZ':
    1834918358        ovar1 = onc.variables['temp']
    … …  
    1835118360            temp = ovar1[:].astype('float64')
    1835218361        else:
    18353             temp = ovar1[:].astype('float')
     18362            temp = ovar1[:]
    1835418363    else:
    1835518364        ovar1 = onc.variables['ta']
    … …  
    1835718366            temp = ovar1[:].astype('float64')
    1835818367        else:
    18359             temp = ovar1[:].astype('float')
    18360 
     18368            temp = ovar1[:]
     18369
     18370    print '  ' + fname + ": Opening output file 'pinterp.nc'"
    1836118371    onewnc = NetCDFFile(ofile, 'w')
    1836218372# Creation of dimensions
    … …  
    1844618456                print '  ' + fname + ': De-staggering x-wind variable !!'
    1844718457                print '    from:', ovarin.shape, 'to', (dimt, dimz, dimy, dimx)
    18448                 varin = np.zeros((dimt, dimz, dimy, dimx), dtype=np.float64)
    18449                 varin[:] = 0.5*(ovarin[:,:,:,0:dimx] + ovarin[:,:,:,1:dimx+1])
     18458                if isFR64:
     18459                    varin = np.zeros((dimt, dimz, dimy, dimx), dtype=np.float64)
     18460                else:
     18461                    varin = np.zeros((dimt, dimz, dimy, dimx), dtype=np.float)
     18462                # Not pro, but less memory problems!
     18463                for it in range(dimt):
     18464                    varin[it,:,::] = 0.5*(ovarin[it,:,:,0:dimx] + ovarin[it,:,:,1:dimx+1])
    1845018465                isgeop = False
    1845118466                varattrs = gen.variables_values('ua')
    … …  
    1845918474                print '  ' + fname + ': De-staggering y-wind variable !!'
    1846018475                print '    from:', ovarin.shape, 'to', (dimt, dimz, dimy, dimx)
    18461                 varin = np.zeros((dimt, dimz, dimy, dimx), dtype=np.float64)
     18476                if isFR64:
     18477                    varin = np.zeros((dimt, dimz, dimy, dimx), dtype=np.float64)
     18478                else:
     18479                    varin = np.zeros((dimt, dimz, dimy, dimx), dtype=np.float)
    1846218480                varin[:] = 0.5*(ovarin[:,:,0:dimy,:] + ovarin[:,:,1:dimy+1,:])
    1846318481                isgeop = False
    … …  
    1851318531                tini = itt
    1851418532                tend = itt + dimtfrac
    18515                 varint = varin[tini:tend,:,:,:].transpose()
    18516                 prest = pres[tini:tend,:,:,:].transpose()
    18517                 psfct = psfc[tini:tend,:,:].transpose()
    18518                 hgtt = hgt.transpose()
    18519                 tempt = temp[tini:tend,:,:,:].transpose()
    18520                 qvt = qv[tini:tend,:,:,:].transpose()
    18521 
    18522                 varinterp = fin.module_forinterpolate.interp( data_in=varint,            \
    18523                   pres_field=prest, interp_levels=interplevs, psfc=psfct,                \
    18524                   ter=hgtt, tk=tempt, qv=qvt, linlog=linloginterp,                       \
    18525                   extrapolate=extrap, geopt=isgeop, missing=gen.fillValueF)
    18526                 newvar[tini:tend,:,:,:] = varinterp.transpose()
     18533                if np.prod([tend-tini+1, dimz, dimx, dimy]) > 150*150*39:
     18534                    dimyfrac = 50
     18535                    print warnmsg
     18536                    print '  ' + fname + ': variable to interpolate:', varin.shape,      \
     18537                      'too big!!'
     18538                    print '    p-interpolation will be done by yaxis-slices of 10 grid-points'
     18539                    print range(0,dimyfrac*int(dimy/dimyfrac),dimyfrac)
     18540                    for iyy in range(0,dimyfrac*int(dimy/dimyfrac),dimyfrac):
     18541                        yini = iyy
     18542                        yend = iyy + dimyfrac
     18543                        varint = varin[tini:tend,:,yini:yend,:].transpose()
     18544                        prest = pres[tini:tend,:,yini:yend,:].transpose()
     18545                        psfct = psfc[tini:tend,yini:yend,:].transpose()
     18546                        hgtt = hgt[yini:yend,:].transpose()
     18547                        tempt = temp[tini:tend,:,yini:yend,:].transpose()
     18548                        qvt = qv[tini:tend,:,yini:yend,:].transpose()
     18549                        Nplevs = len(interplevs)
     18550                        Ndt = len(range(tini,tend))
     18551                        Ndy = len(range(yini,yend))
     18552                        print 'Lluis before!!!!'
     18553                        varinterp = fin.module_forinterpolate.interp( data_in=varint,\
     18554                          pres_field=prest, interp_levels=interplevs, psfc=psfct,    \
     18555                          ter=hgtt, tk=tempt, qv=qvt, linlog=linloginterp,           \
     18556                          extrapolate=extrap, geopt=isgeop, missing=gen.fillValueF,  \
     18557                          ix=dimx, iy=Ndy, iz=dimz, it=Ndt, num_metgrid_levels=Nplevs)
     18558                        print 'Lluis done!!!!'
     18559                        newvar[tini:tend,:,yini:yend,:] = varinterp.transpose()
     18560                else:
     18561                    varint = varin[tini:tend,:,:,:].transpose()
     18562                    prest = pres[tini:tend,:,:,:].transpose()
     18563                    psfct = psfc[tini:tend,:,:].transpose()
     18564                    hgtt = hgt.transpose()
     18565                    tempt = temp[tini:tend,:,:,:].transpose()
     18566                    qvt = qv[tini:tend,:,:,:].transpose()
     18567                    Nplevs = len(interplevs)
     18568                    Ndt = len(range(tini,tend))
     18569                    print 'Lluis before!!!!'
     18570                    varinterp = fin.module_forinterpolate.interp( data_in=varint,    \
     18571                      pres_field=prest, interp_levels=interplevs, psfc=psfct,        \
     18572                      ter=hgtt, tk=tempt, qv=qvt, linlog=linloginterp,               \
     18573                      extrapolate=extrap, geopt=isgeop, missing=gen.fillValueF,      \
     18574                      ix=dimx, iy=dimy, iz=dimz, it=Ndt, num_metgrid_levels=Nplevs)
     18575                    print 'Lluis done!!!!'
     18576                    newvar[tini:tend,:,:,:] = varinterp.transpose()
    1852718577
    1852818578            if dimtfrac != dimt and tend != dimt-1:
Note: See TracChangeset for help on using the changeset viewer.