Changeset 2372 in lmdz_wrf

Timestamp:

Feb 26, 2019, 9:54:13 PM (6 years ago)

Author:

lfita

Message:

Adding new version of 'pinterp' using `range_slicing' !!

Location:

trunk/tools

Files:

• generic_tools.py (modified) (2 diffs)
• nc_var_tools.py (modified) (1 diff)
• variables_values.dat (modified) (2 diffs)

trunk/tools/generic_tools.py

@@ -15227 +15227 @@
 
     return loops1D
+#print Nloops_1D(['t', 'z', 'y', 'x'], [56, 39, 99, 99], {'y': 99, 't': 56})
+#quit()
 
 def range_slicing(dimns, dimvs, intdims):
@@ -15309 +15311 @@
 
     # Getting all the combinations of the interval slices
-    intcombs = Nloops_1D(intloopvals)
-    Tslicesize = intcombs.shape[0]
+    if np.all(np.array(intloopvals) == 0):
+        print warnmsg
+        print '  ' + fname + ': not actual intervals for any dimension !!'
+        for dimn in intdims.keys():
+            print '  intervals ', dimn, ':', intslices[dimn] 
+        print '    creation of empty intervals ...'
+        Nintdims = len(list(intdims.keys()))
+        intcombs = [list(np.zeros((Nintdims), dtype=int))]
+        for dimn in dimns:
+            if intdims.has_key(dimn): 
+                intslices[dimn] = [0, dims[dimn]]
+        Tslicesize = 1
+    else:
+        intcombs = Nloops_1D(intloopvals)
+        Tslicesize = intcombs.shape[0]
 
     #print '  ' + fname + ': Total number of slices to provide:', Tslicesize

trunk/tools/nc_var_tools.py

@@ -18809 +18809 @@
     print '  ' + fname + ' interpolating at pressure levels _______' 
     print '     ', interplevs
+    usaid = False
+    uersaid = False
+    vsaid = False
+    versaid = False
+    tsaid = False
+    ysaid = False
+    wsaid = False
+
+    Nplevs = len(interplevs)
+    for vn in varns:
+        print "      '" + vn + "' ..."
+        newvarattr = {}
+        varin = None
+
+        # FORCING the splitting
+#        # There is a memroy issue, thus is necessary to split the matrix...
+#        if np.prod(list(temp.shape))*8 > 3*gen.oneGB:
+#            #dimtfrac = np.min([5,dimt])
+#            if not tsaid:
+#                print warnmsg
+#                print '  ' + fname + ': domain to interpolate:', temp.shape,         \
+#                  'too big!!'
+#                print '   p-interpolation will be done by time-slices of', dimtfrac, \
+#                  'time-steps'
+#                print range(0,dimtfrac*int(dimt/dimtfrac),dimtfrac)
+#        else:
+#            dimtfrac = dimt
+#        if np.prod([dimz, dimx, dimy])*8 > gen.oneGB:
+#            #dimyfrac = 50
+#            if not ysaid:
+#                print warnmsg
+#                print '  ' + fname + ': variable to interpolate:', temp.shape,       \
+#                  'too big!!'
+#                print '    p-interpolation will be done by yaxis-slices of',         \
+#                  dimyfrac,' grid-points'
+#                print range(0,dimyfrac*int(dimy/dimyfrac),dimyfrac)
+#        else:
+#            dimyfrac = dimy
+        if not tsaid:
+            print warnmsg
+            print '   p-interpolation will be done by time-slices of', dimtfrac,     \
+              'time-steps'
+            print range(0,dimtfrac*int(dimt/dimtfrac),dimtfrac)
+        if not ysaid:
+            print warnmsg
+            print '    p-interpolation will be done by yaxis-slices of',             \
+              dimyfrac,' grid-points'
+            print range(0,dimyfrac*int(dimy/dimyfrac),dimyfrac)
+
+        # Getting slices
+        tzslices, tzinds = gen.range_slicing(['t', 'z', 'y', 'x'], list(temp.shape), \
+          {'t': dimtfrac, 'y': dimyfrac})
+
+        varinterp = np.zeros(tuple([dimx,dimyfrac,len(interplevs),dimtfrac]), dtype=np.float)
+        TOTslices = len(tzslices)
+
+        # Splitting
+        for islc in range(TOTslices):
+            slcv = tzslices[islc]
+            slci = tzinds[islc]
+            tini = slci[0,0]
+            tend = slci[0,1]
+            yini = slci[2,0]
+            yend = slci[2,1]
+
+            Ndy = len(range(yini,yend))
+            Ndt = len(range(tini,tend))
+            #Slice tuple
+            tslc = tuple([slice(tini,tend), slice(0,dimz), slice(yini,yend),         \
+              slice(0,dimx)])
+            tslcx1 = tuple([slice(tini,tend), slice(0,dimz), slice(yini,yend),       \
+              slice(1,dimx+1)])
+            tslcy1= tuple([slice(tini,tend), slice(0,dimz), slice(yini+1,yend+1),    \
+              slice(0,dimx)])
+
+            if gen.searchInlist(notCHK, vn):
+                if vn == 'WRFght': 
+                    varin = geop[tslc]
+                    isgeop = True
+                    varattrs = gen.variables_values('WRFght')
+                    CFvn = varattrs[0]
+                    newvarattr['standard_name'] = varattrs[1]
+                    newvarattr['long_name'] = varattrs[4].replace('|',' ')
+                    newvarattr['units'] = varattrs[5]
+                elif vn == 'WRFhus': 
+                    # specific humidity
+                    varin = qv[tslc]/(1.+qv[tslc])
+                    isgeop = False
+                    varattrs = gen.variables_values('WRFhus')
+                    CFvn = varattrs[0]
+                    newvarattr['standard_name'] = varattrs[1]
+                    newvarattr['long_name'] = varattrs[4].replace('|',' ')
+                    newvarattr['units'] = varattrs[5]
+                elif vn == 'WRFt':
+                    varin = temp[tslc]
+                    isgeop = False
+                    varattrs = gen.variables_values('WRFt')
+                    CFvn = varattrs[0]
+                    newvarattr['standard_name'] = varattrs[1]
+                    newvarattr['long_name'] = varattrs[4].replace('|',' ')
+                    newvarattr['units'] = varattrs[5]
+                elif vn == 'WRFu':
+                    ovarin = onc.variables['U']
+                    if not usaid:
+                        print infmsg
+                        print '  ' + fname + ': De-staggering x-wind variable !!'
+                        print '    from:',ovarin[tini:tend,:,yini:yend,0:dimx].shape,\
+                          'to', (Ndt,dimz,Ndy,dimx)
+                        usaid = True
+                    if isFR64:
+                        varin = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float64)
+                    else:
+                        varin = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float)
+                    varin[:] = 0.5*(ovarin[tslc] + ovarin[tslcx1])
+                    isgeop = False
+                    varattrs = gen.variables_values('ua')
+                    CFvn = varattrs[0]
+                    newvarattr['standard_name'] = varattrs[1]
+                    newvarattr['long_name'] = varattrs[4].replace('|',' ')
+                    newvarattr['units'] = varattrs[5]
+                elif vn == 'WRFuer':
+                    ovarin = onc.variables['U']
+                    ovarin2 = onc.variables['V']
+                    # Earth-rotating
+                    osina =  onc.variables['SINALPHA']
+                    ocosa =  onc.variables['COSALPHA']
+                    if not uersaid:
+                        print infmsg
+                        print '  ' + fname + ': De-staggering x-wind variable & ' +  \
+                          'Earth rotating!!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        uersaid = True
+                    if isFR64:
+                        varin = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float64)
+                        varin0 = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float64)
+                        varin02 = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float64)
+                    else:
+                        varin = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float)
+                        varin0 = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float)
+                        varin02 = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float)
+                    for it in range(tini,tend):
+                        varin0[it-tini,:,:,:]=0.5*(ovarin[it,:,yini:yend,0:dimx] +   \
+                          ovarin[it,:,yini:yend,1:dimx+1])
+                        varin02[it-tini,:,:,:] = 0.5*(ovarin2[it,:,yini:yend,:] +    \
+                          ovarin2[it,:,yini+1:yend+1,:])
+                    for iz in range(dimz):
+                        varin[:,iz,:,:] = varin0[:,iz,:,:]*                          \
+                          ocosa[tini:tend,yini:yend,:] - varin02[:,iz,:,:]*          \
+                          osina[tini:tend,yini:yend,:]
+                    isgeop = False
+                    varattrs = gen.variables_values('uer')
+                    CFvn = varattrs[0]
+                    newvarattr['standard_name'] = varattrs[1]
+                    newvarattr['long_name'] = varattrs[4].replace('|',' ')
+                    newvarattr['units'] = varattrs[5]
+                elif vn == 'WRFv':
+                    ovarin = onc.variables['V']
+                    if not vsaid:
+                        print infmsg
+                        print '  ' + fname + ': De-staggering y-wind variable !!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        vsaid = True
+                    if isFR64:
+                        varin = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float64)
+                    else:
+                        varin = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float)
+                    # Not pro, but less memory problems!
+                    for it in range(tini,tend):
+                        varin[it-tini,:,:,:] = 0.5*(ovarin[it,:,yini:yend,:] +       \
+                          ovarin[it,:,yini+1:yend+1,:])
+                    isgeop = False
+                    varattrs = gen.variables_values('va')
+                    CFvn = varattrs[0]
+                    newvarattr['standard_name'] = varattrs[1]
+                    newvarattr['long_name'] = varattrs[4].replace('|',' ')
+                    newvarattr['units'] = varattrs[5]
+                elif vn == 'WRFver':
+                    ovarin = onc.variables['U']
+                    ovarin2 = onc.variables['V']
+                    # Earth-rotating
+                    osina =  onc.variables['SINALPHA']
+                    ocosa =  onc.variables['COSALPHA']
+                    if not versaid:
+                        print infmsg
+                        print '  ' + fname + ': De-staggering y-wind variable & ' +  \
+                          'Earth rotating!!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        versaid = True
+                    if isFR64:
+                        varin = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float64)
+                        varin0 = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float64)
+                        varin02 = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float64)
+                    else:
+                        varin = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float)
+                        varin0 = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float)
+                        varin02 = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float)
+                    # Not pro, but less memory problems!
+                    for it in range(tini,tend):
+                        varin0[it-tini,:,:,:]=0.5*(ovarin[it,:,yini:yend,0:dimx] +   \
+                          ovarin[it,:,yini:yend,1:dimx+1])
+                        varin02[it-tini,:,:,:] = 0.5*(ovarin2[it,:,yini:yend,:] +    \
+                           ovarin2[it,:,yini+1:yend+1,:])
+                    for iz in range(dimz):
+                        varin[:,iz,:,:] = varin0[:,iz,:,:]*                          \
+                          osina[tini:tend,yini:yend,:] + varin02[:,iz,:,:]*          \
+                          ocosa[tini:tend,yini:yend,:]
+                    isgeop = False
+                    varattrs = gen.variables_values('ver')
+                    CFvn = varattrs[0]
+                    newvarattr['standard_name'] = varattrs[1]
+                    newvarattr['long_name'] = varattrs[4].replace('|',' ')
+                    newvarattr['units'] = varattrs[5]
+                elif vn == 'WRFw':
+                    ovarin = onc.variables['W']
+                    if not wsaid:
+                        print infmsg
+                        print '  ' + fname + ': De-staggering z-wind variable !!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        wsaid = True
+                    if isFR64:
+                        varin = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float64)
+                    else:
+                        varin = np.zeros((Ndt,dimz,Ndy,dimx), dtype=np.float)
+                    # Not pro, but less memory problems!
+                    for it in range(tini,tend):
+                        varin[it-tini,:,:,:]=0.5*(ovarin[it,0:dimz,yini:yend,:] +    \
+                          ovarin[it,1:dimz+1,yini:yend,:])
+                    isgeop = False
+                    varattrs = gen.variables_values('wa')
+                    CFvn = varattrs[0]
+                    newvarattr['standard_name'] = varattrs[1]
+                    newvarattr['long_name'] = varattrs[4].replace('|',' ')
+                    newvarattr['units'] = varattrs[5]
+                elif vn == 'pres':
+                    varin = pres[tslc]
+                    isgeop = False
+                    varattrs = gen.variables_values('pres')
+                    CFvn = varattrs[0]
+                    newvarattr['standard_name'] = varattrs[1]
+                    newvarattr['long_name'] = varattrs[4].replace('|',' ')
+                    newvarattr['units'] = varattrs[5]
+#            elif not gen.searchInlist(MODvarrequired, vn) or vn == 'QVAPOR':
+            else:
+                ovarin = onc.variables[vn]
+                varinf = variable_inf(ovarin)
+                varin = ovarin[tslc]
+                # de-staggering
+                #print '  Lluis:', ovarin.dimensions,'unsDIM',  unstaggerDIM, 'CFdim:', CFdimvs
+                #if gen.searchInlist(ovarin.dimensions, unstaggerDIM):
+                #    varin = gen.stagger_unstagger(varin,list(ovarin.dimensions),MODdimvs)
+                #else:
+                #    varin = gen.stagger_unstagger(varin,list(ovarin.dimensions),CFdimvs)
+                isgeop = False
+                if vn == 'z': isgeop = True
+                varattrs = gen.variables_values(vn)
+                CFvn = varattrs[0]
+                for attrn in ovarin.ncattrs():
+                    attrv = ovarin.getncattr(attrn)
+                    newvarattr[attrn] = attrv
+                newvarattr['standard_name'] = varattrs[1]
+                newvarattr['long_name'] = varattrs[4].replace('|',' ')
+
+            if varin is not None:
+                if not gen.searchInlist(onewnc.variables,CFvn): 
+                    newvar = onewnc.createVariable(CFvn, 'f4', tuple(newMODdims),    \
+                      fill_value=gen.fillValueF)
+                varint = varin.transpose()
+                prest = pres[tslc].transpose()
+                psfct = psfc[tini:tend,yini:yend,:].transpose()
+                hgtt = hgt[yini:yend,:].transpose()
+                tempt = temp[tslc].transpose()
+                qvt = qv[tslc].transpose()
+                varinterp = fin.module_forinterpolate.interp( data_in=varint,        \
+                  pres_field=prest, interp_levels=interplevs, psfc=psfct,            \
+                  ter=hgtt, tk=tempt, qv=qvt, linlog=linloginterp,                   \
+                  extrapolate=extrap, geopt=isgeop, missing=gen.fillValueF,          \
+                  ix=dimx, iy=Ndy, iz=dimz, it=Ndt, num_metgrid_levels=Nplevs)
+                newvar[tini:tend,:,yini:yend,:] = varinterp.transpose()
+                onewnc.sync()
+
+        for attrn in newvarattr.keys():
+            if attrn != '_FillValue':
+                attr = set_attribute(newvar, attrn, newvarattr[attrn])
+        onewnc.sync()
+
+    # Global attributes
+    add_global_PyNCplot(onewnc, 'nc_var_tools.py', fname, '2.0')
+
+    onewnc.setncattr('oringal', 'Interpolation using p_interp.F90 subroutines')
+    newattr = set_attributek(onewnc, 'original_subroutines_author', unicode('Cindy '+\
+      'Bruy' + unichr(232) + 're'), 'U')
+    onewnc.setncattr('original_subroutines_date', 'November 2007')
+
+    for attrn in onc.ncattrs():
+        attrv = onc.getncattr(attrn)
+        attr = set_attribute(onewnc, attrn, attrv)    
+
+    onc.close()
+
+    onewnc.sync()
+    onewnc.close()
+
+    print fname + ": successfull written of '" + ofile + "' !!"
+
+    return
+
+
+def pinterp_old(values, ncfile, variables):
+    """ Function to vertically interpolate using subroutines from the p_interp.F90 NCAR program
+      Using fortran codes: module_generic.F90
+      values= [interplevs],[linloginterp],[extrap],[dimtfrac],[dimyfrac]
+        [interplevs]: ':' separated list of pressure values to interpolate [Pa] (decreassing in Pa)
+        [linloginterp]: kind of linear interpolation
+          1: linear interp in pressure
+          2: linear interp in ln(pressure)
+        [extrap]: whether to set to missing value below/above model ground and top (0), or extrapolate (1)
+        [dimtfrac]: fractions of time to interpolate to avoid memory issues ('auto' for 5)
+        [dimyfrac]: fractions of y-axis to interpolate to avoid memory issues ('auto' for 50)
+      ncfile= WRF file to use
+      variables = ',' separated list of names of variables to interpolate ('all', fo all 4D-atmospheric variables)
+        'WRFght': for WRF geopotential height
+        'WRFrh': for WRF relative humidity
+        'WRFt': for WRF temperature
+        'WRFu': for WRF x-wind de-staggered
+        'WRFuer': for WRF x-wind de-staggered Earth-rotated
+        'WRFv': for WRF y-wind de-staggered
+        'WRFver': for WRF y-wind de-staggered Earth-rotated
+        'WRFw': for WRF z-wind de-staggered
+    """
+    import module_ForInt as fin
+    fname = 'pinterp'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print pinterp.__doc__
+        quit()
+
+    arguments = '[interplevs],[linloginterp],[extrap],[dimtfrac],[dimyfrac]'
+    gen.check_arguments(fname, values, arguments, ',')
+
+    if isFR64:
+        interplevs = np.array(values.split(',')[0].split(':'), dtype=np.float64)
+    else:
+        interplevs = np.array(values.split(',')[0].split(':'), dtype=np.float)
+    linloginterp = np.int32(values.split(',')[1])
+    extrap = np.int32(values.split(',')[2])
+    dimtfrac = gen.auto_val(values.split(',')[3], 5)
+    dimyfrac = gen.auto_val(values.split(',')[4], 50)
+
+    ofile = 'pinterp.nc'
+
+    CFdims = ['time', 'pres', 'lat', 'lon']
+    # Just in case there are weird units, attribute value will be taken from input file
+    CFattrs = ['standard_name', 'long_name']
+
+    onc = NetCDFFile(ncfile, 'r')
+
+    # Gessing orgin of the file
+    dimsinfile = onc.dimensions.keys()
+    varsinfile = onc.variables.keys()
+ 
+    if gen.searchInlist(dimsinfile,'bottom_top'):
+        modname = 'WRF'
+        print warnmsg
+        print '  ' + fname + ": gessing that file '" + ncfile + "' comes from '" +   \
+          modname + "' !!"
+
+        # Variables to interpolate
+        MODdims = ['Time', 'bottom_top', 'south_north', 'west_east']
+        newMODdims = ['Time', 'pres', 'south_north', 'west_east']
+        notCHK = ['WRFght', 'WRFhus', 'WRFrh', 'WRFt', 'WRFu', 'WRFuer', 'WRFv',     \
+          'WRFver', 'WRFw']
+        MODvarrequired = ['P', 'PB', 'PSFC', 'PH', 'PHB', 'HGT', 'T', 'QVAPOR',      \
+          'XLONG', 'XLAT', 'Times']
+        MODvardims = ['XLONG', 'XLAT', 'Times']
+    elif gen.searchInlist(dimsinfile,'presnivs'):
+        modname = 'LMDZ'
+        print warnmsg
+        print '  ' + fname + ": gessing that file '" + ncfile + "' comes from '" +   \
+          modname + "' !!"
+
+        # Variables to interpolate
+        MODdims = ['time_counter', 'presnivs', 'lat', 'lon']
+        newMODdims = ['time_counter', 'pres', 'lat', 'lon']
+        notCHK = []
+        MODvarrequired = ['pres', 'psol', 'geop', 'phis', 'temp', 'ovap', 'lon',     \
+          'lat', 'time_counter']
+        MODvardims = ['lon', 'lat', 'time_counter']
+    elif gen.searchInlist(dimsinfile,'pres') and gen.searchInlist(varsinfile,        \
+      'psol'):
+        modname = 'cfLMDZ'
+        print warnmsg
+        print '  ' + fname + ": gessing that file '" + ncfile + "' comes from '" +   \
+          modname + "' ('pres, time_counter' variables renamed as 'p, time') !!"
+
+        # Variables to interpolate
+        MODdims = ['time', 'pres', 'lat', 'lon']
+        newMODdims = ['time', 'pres', 'lat', 'lon']
+        notCHK = []
+        MODvarrequired = ['p', 'psol', 'geop', 'phis', 'temp', 'ovap', 'lon',     \
+          'lat', 'time']
+        MODvardims = ['lon', 'lat', 'time']
+    else:
+        modname = 'CF'
+        print warnmsg
+        print '  ' + fname + ": gessing that file '" + ncfile + "' comes from '" +   \
+          modname + "' !!"
+
+        # Variables to interpolate
+        MODdims = ['time', 'pres', 'lat', 'lon']
+        newMODdims = ['time', 'pres', 'lat', 'lon']
+        notCHK = []
+        MODvarrequired = ['p', 'ps', 'z', 'orog', 'ta', 'hus', 'lon', 'lat', 'time']
+        MODvardims = ['lon', 'lat', 'time']
+
+    varns = gen.str_list(variables, ',')
+    if variables == 'all':
+        varns = []
+        for vn in onc.variables:
+            ovn = onc.variables[vn]
+            coincdims = set(WRFdims) & set(ovn.dimensions)
+           
+            if len(MODdims) == len(ovn.dimensions) and len(coincdims) == len(MODdims):
+                varns.append(vn)
+        varns = varns + notCHK
+        
+    for vn in varns:
+        if not gen.searchInlist(onc.variables, vn) and                               \
+          not gen.searchInlist(notCHK, vn):
+            print errormsg
+            print '  ' + fname + "': file '" + ncfile + "' does not have variable '"+\
+              vn + "' !!"
+            quit(-1)
+
+    # looking for model required variables
+    for var in MODvarrequired:
+        if not gen.searchInlist(onc.variables, var):
+            print errormsg
+            print '  ' + fname + ": file '" + ncfile + "' does not have required " + \
+              "variable '" + var + "' !!"
+            quit(-1)
+
+# Computing necessary variables
+    # pressure
+    if modname == 'WRF':
+        ovar1 = onc.variables['P']
+        ovar2 = onc.variables['PB']
+
+        pres0 = ovar1[:] + ovar2[:]
+        if isFR64:
+            pres = pres0.astype('float64')
+        else:
+            pres = pres0.copy()
+
+        dimx = pres.shape[3]
+        dimy = pres.shape[2]
+        dimz = pres.shape[1]
+        dimt = pres.shape[0]
+        MODdimvs = {'Time':dimt, 'bottom_top':dimz, 'south_north':dimy,              \
+          'west_east': dimx}
+        CFdimvs = {'time': dimt, 'bottom_top': dimz, 'lat': dimy, 'lon': dimx}
+        unstaggerDIM = 'west_east'
+    elif modname == 'LMDZ':
+        ovar1 = onc.variables['pres']
+        if isFR64:
+            pres = ovar1[:].astype('float64')
+        else:
+            pres = ovar1[:]
+        dimx = pres.shape[3]
+        dimy = pres.shape[2]
+        dimz = pres.shape[1]
+        dimt = pres.shape[0]
+        MODdimvs = {'time_counter':dimt, 'presnivs':dimz, 'lat':dimy, 'lon': dimx}
+        CFdimvs = {'time': dimt, 'presnivs': dimz, 'lat': dimy, 'lon': dimx}
+        unstaggerDIM = 'lon'
+    elif modname == 'cfLMDZ':
+        ovar1 = onc.variables['p']
+        if isFR64:
+            pres = ovar1[:].astype('float64')
+        else:
+            pres = ovar1[:]
+        dimx = pres.shape[3]
+        dimy = pres.shape[2]
+        dimz = pres.shape[1]
+        dimt = pres.shape[0]
+        MODdimvs = {'time':dimt, 'pres':dimz, 'lat':dimy, 'lon': dimx}
+        CFdimvs = {'time': dimt, 'pres': dimz, 'lat': dimy, 'lon': dimx}
+        unstaggerDIM = 'lon'
+    else:
+        ovar1 = onc.variables['p']
+        if isFR64:
+            pres = ovar1[:].astype('float64')
+        else:
+            pres = ovar1[:]
+        dimx = pres.shape[3]
+        dimy = pres.shape[2]
+        dimz = pres.shape[1]
+        dimt = pres.shape[0]
+        MODdimvs = {'time':dimt, 'pres':dimz, 'lat':dimy, 'lon': dimx}
+        CFdimvs = {'time': dimt, 'pres': dimz, 'lat': dimy, 'lon': dimx}
+        unstaggerDIM = 'lon'
+
+    print '  ' + fname + ': CF dimension lengths:', CFdimvs
+
+    # sfc pressure
+    if modname == 'WRF':
+        ovar1 = onc.variables['PSFC']
+        if isFR64:
+            psfc = ovar1[:].astype('float64')
+        else:
+            psfc = ovar1[:]
+    elif modname == 'LMDZ' or modname == 'cfLMDZ':
+        ovar1 = onc.variables['psol']
+        if isFR64:
+            psfc = ovar1[:].astype('float64')
+        else:
+            psfc = ovar1[:]
+    else:
+        ovar1 = onc.variables['ps']
+        if isFR64:
+            psfc = ovar1[:].astype('float64')
+        else:
+            psfc = ovar1[:]
+
+    # geopotential height
+    if modname == 'WRF':
+        ovar1 = onc.variables['PH']
+        ovar2 = onc.variables['PHB']
+        geop0 = ovar1[:] + ovar2[:]
+        unstg = list(geop0.shape)
+        unstg[1] = unstg[1] - 1
+        if isFR64:
+            geop = np.zeros(tuple(unstg), dtype=np.float64)
+            geop = 0.5*(geop0[:,1:dimz+1,:,:] + geop0[:,0:dimz,:,:]).astype('float64')
+        else:
+            geop = np.zeros(tuple(unstg), dtype=np.float)
+            geop = 0.5*(geop0[:,1:dimz+1,:,:] + geop0[:,0:dimz,:,:])
+    elif modname == 'LMDZ' or modname == 'cfLMDZ':
+        ovar1 = onc.variables['geop']
+        if isFR64:
+            geop = ovar1[:].astype('float64')
+        else:
+            geop = ovar1[:]
+    else:
+        ovar1 = onc.variables['z']
+        if isFR64:
+            geop = ovar1[:].astype('float64')
+        else:
+            geop = ovar1[:]
+
+    # terrain height
+    if modname == 'WRF':
+        ovar1 = onc.variables['HGT']
+        if isFR64:
+            hgt = ovar1[0,:,:].astype('float64')
+        else:
+            hgt = ovar1[0,:,:]
+    elif modname == 'LMDZ' or modname == 'cfLMDZ':
+        grav = 9.81
+        ovar1 = onc.variables['phis']
+        if isFR64:
+            hgt = (ovar1[0,:,:]/grav).astype('float64')
+        else:
+            hgt = (ovar1[0,:,:]/grav)
+    else:
+        ovar1 = onc.variables['orog']
+        if isFR64:
+            hgt = ovar1[:].astype('float64')
+        else:
+            hgt = ovar1[:]
+
+    # water vapour mixing ratio
+    if modname == 'WRF':
+        ovar1 = onc.variables['QVAPOR']
+        if isFR64:
+            qv = ovar1[:].astype('float64')
+        else:
+            qv = ovar1[:]
+    elif modname == 'LMDZ' or modname == 'cfLMDZ':
+        ovar1 = onc.variables['ovap']
+        if isFR64:
+            qv = ovar1[:].astype('float64')
+        else:
+            qv = ovar1[:]
+    else:
+        ovar1 = onc.variables['hus']
+        if isFR64:
+            qv = ovar1[:].astype('float64')
+        else:
+            qv = ovar1[:]
+
+    # temperature
+    if modname == 'WRF':
+        if isFR64:
+            Rd = np.float64(287.04)
+            Cp = np.float64(7.*Rd/2.)
+            RCP = np.float64(Rd/Cp)
+            p0 = np.float64(100000.)
+        else:
+            Rd = 287.04
+            Cp = 7.*Rd/2.
+            RCP = Rd/Cp
+            p0 = 100000.
+        ovar10 = onc.variables['T']
+        var10 = ovar10[:]
+        if isFR64:
+            ovar1 = (var10).astype('float64')
+        else:
+            ovar1 = var10[:]
+        temp0 = (ovar1[:]+300.)*(pres[:]/p0)**RCP
+        if isFR64:
+            temp = temp0.astype('float64')
+        else:
+            temp = temp0.copy()
+    elif modname == 'LMDZ' or modname == 'cfLMDZ':
+        ovar1 = onc.variables['temp']
+        if isFR64:
+            temp = ovar1[:].astype('float64')
+        else:
+            temp = ovar1[:]
+    else:
+        ovar1 = onc.variables['ta']
+        if isFR64:
+            temp = ovar1[:].astype('float64')
+        else:
+            temp = ovar1[:]
+
+    print '  ' + fname + ": Opening output file 'pinterp.nc'"
+    onewnc = NetCDFFile(ofile, 'w')
+# Creation of dimensions
+    if modname == 'WRF':
+        newdim = onewnc.createDimension('west_east', dimx)
+        newdim = onewnc.createDimension('south_north', dimy)
+        newdim = onewnc.createDimension('pres', len(interplevs))
+        newdim = onewnc.createDimension('Time', None)
+    elif modname == 'LMDZ':
+        newdim = onewnc.createDimension('lon', dimx)
+        newdim = onewnc.createDimension('lat', dimy)
+        newdim = onewnc.createDimension('pres', len(interplevs))
+        newdim = onewnc.createDimension('time_counter', None)
+    else:
+        newdim = onewnc.createDimension('lon', dimx)
+        newdim = onewnc.createDimension('lat', dimy)
+        newdim = onewnc.createDimension('pres', len(interplevs))
+        newdim = onewnc.createDimension('time', None)
+
+# Creation of variable dimensions
+    for var in MODvardims:
+        ovar = onc.variables[var]
+        varinf = variable_inf(ovar)
+        newvard = []
+        for vd in varinf.dimns:
+            if not gen.searchInlist(onewnc.dimensions, vd) and                       \
+              not gen.searchInlist(CFdims, vd):
+                newdim = onewnc.createDimension(vd, len(onc.dimensions[vd]))
+            if gen.searchInlist(CFdims,vd):
+                icfdim = CFdims.index(vd)
+                newvard.append(MODdims[icfdim])
+            else:
+                newvard.append(vd)
+
+        newvar = onewnc.createVariable(var, nctype(varinf.dtype),                    \
+          tuple(newvard))
+        newvar[:] = ovar[:]
+        for attrn in ovar.ncattrs():
+            attrv = ovar.getncattr(attrn)
+            attr = set_attribute(newvar, attrn, attrv)    
+
+    # Creation of pressure variable dimension
+    newvar = onewnc.createVariable('pres', 'f8', ('pres'))
+    newvar[:] = interplevs[:]
+    basicvardef(newvar, 'pressure', 'Pressure', 'Pa')
+    attr = set_attribute(newvar, 'positive', 'down')
+    attr = set_attribute(newvar, 'axis', 'Z')
+
+    print '  ' + fname + ' interpolating at pressure levels _______' 
+    print '     ', interplevs
     tsaid = False
     ysaid = False

trunk/tools/variables_values.dat

@@ -545 +545 @@
 Zonal wind, ua, eastward_wind, -30., 30., eastward|wind, ms-1, seismic, $ua$, ua
 LVITU, ua, eastward_wind, -30., 30., eastward|wind, ms-1, seismic, $ua$, ua
+uer, uer, eastward_wind_er, -30., 30., Earth|rotated|eastward|wind, ms-1, seismic, $uer$, uer
 uas, uas, eastward_wind, -30., 30., eastward|10m|wind, ms-1, seismic, $uas$, uas
 UAS, uas, eastward_wind, -30., 30., eastward|10m|wind, ms-1, seismic, $uas$, uas
@@ -558 +559 @@
 Meridional wind, va, northward_wind, -30., 30., northward|wind, ms-1, seismic, $va$, va
 LVITV, va, northward_wind, -30., 30., northward|wind, ms-1, seismic, $va$, va
+ver, ver, northward_wind_er, -30., 30., Earth|rotated|northward|wind, ms-1, seismic, $ver$, ver
 valley, valley, valley_point, 0., 1., valley|grd|point, 1, Binary, $valley$, valley
 vas, vas, northward_wind, -30., 30., northward|10m|wind, ms-1, seismic, $vas$, vas