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

Timestamp:

Mar 20, 2018, 8:22:37 PM (7 years ago)

Author:

lfita

Message:

Adding a new 'pinterp' with

• t/y-splitting on variable input deifinition
• dimtfrac and dimyfrac now as paramters

File:

• trunk/tools/nc_var_tools.py (modified) (2 diffs)

trunk/tools/nc_var_tools.py

@@ -18060 +18060 @@
     """ 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
+
+    Nplevs = len(interplevs)
+    for vn in varns:
+        print "      '" + vn + "' ..."
+        newvarattr = {}
+        varin = None
+
+        # There is a memroy issue, thus is necessary to split the matrix...
+        if np.prod(list(temp.shape)) > 24*150*150*39:
+            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 5 time-steps'
+                print range(0,dimtfrac*int(dimt/dimtfrac),dimtfrac)
+        else:
+            dimtfrac = dimt
+        if np.prod([dimz, dimx, dimy]) > 150*150*39:
+            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
+
+        varinterp = np.zeros(tuple([dimx,dimyfrac,len(interplevs),dimtfrac]), dtype=np.float)
+
+        # Splitting
+        for itt in range(0,dimtfrac*int(dimt/dimtfrac),dimtfrac):
+            tini = itt
+            tend = itt + dimtfrac
+            Ndt = len(range(tini,tend))
+            for iyy in range(0,dimyfrac*int(dimy/dimyfrac),dimyfrac):
+                yini = iyy
+                yend = iyy + dimyfrac
+                Ndy = len(range(yini,yend))
+
+                if gen.searchInlist(notCHK, vn):
+                    if vn == 'WRFght': 
+                        varin = geop[tini:tend,:,yini:yend,:]
+                        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[tini:tend,:,yini:yend,:]/(1.+qv[tini:tend,:,yini:yend,:])
+                        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[tini:tend,:,yini:yend,:]
+                        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']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering x-wind variable !!'
+                        print '    from:', ovarin.shape, 'to', (dimt, dimz, dimy, dimx)
+                        if isFR64:
+                            varin = np.zeros((Ndt, dimz, Ndy, dimx), dtype=np.float64)
+                        else:
+                            varin = np.zeros((Ndt, dimz, Ndy, dimx), dtype=np.float)
+                        varin[tini:tend,:,yini:yend,:] =                             \
+                          0.5*(ovarin[tini:tend,:,yini:yend,0:dimx] +                \
+                          ovarin[tini:tend,:,yini:yend,1:dimx+1])
+
+                        # Not pro, but less memory problems!
+                        #for it in range(dimt):
+                        #    varin[it,:,::] = 0.5*(ovarin[it,:,:,0:dimx] + ovarin[it,:,:,1:dimx+1])
+                        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']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering x-wind variable & ' +  \
+                          'Earth rotating!!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        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,:,yini:yend,:] =                               \
+                              0.5*(ovarin[it,:,yini:yend,0:dimx] +                   \
+                              ovarin[it,:,yini:yend,1:dimx+1])
+                            varin02[it,:,yini:yend,:] =                              \
+                              0.5*(ovarin2[it,:,yini:yend,:] +                       \
+                              ovarin2[it,:,yini+1:yend+1,:])
+                        for iz in range(dimz):
+                            varin[tini:tend,iz,yini:yend,:] =                        \
+                              varin0[tini:tend,iz,yini:yend,:]*                      \
+                              ocosa[tini:tend,yini:yend,:] -                         \
+                              varin02[tini:tend,iz,yini:yend,:]*                     \
+                              osina[tini:tend,yini:yend,:]
+                        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 == 'WRFv':
+                        ovarin = onc.variables['V']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering y-wind variable !!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        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,:,yini:yend,:] = 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']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering y-wind variable & ' +  \
+                          'Earth rotating!!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        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,:,yini:yend,:] =                               \
+                              0.5*(ovarin[it,:,yini:yend,0:dimx] +                   \
+                              ovarin[it,:,yini:yend,1:dimx+1])
+                            varin02[it,:,yini:yend,:] =                              \
+                              0.5*(ovarin2[it,:,yini:yend,:] +                       \
+                              ovarin2[it,:,yini+1:yend+1,:])
+                        for iz in range(dimz):
+                            varin[tini:tend,iz,yini:yend,:] =                        \
+                              varin0[tini:tend,iz,yini:yend,:]*                      \
+                              osina[tini:tend,yini:yend,:] +                         \
+                              varin02[tini:tend,iz,yini:yend,:]*                     \
+                              ocosa[tini:tend,yini:yend,:]
+                        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 == 'WRFw':
+                        ovarin = onc.variables['W']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering z-wind variable !!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        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,:,yini:yend,:] =                                \
+                              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[tini:tend,:,yini:yend,:]
+                        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[tini:tend,:,yini:yend,:]
+                    # 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[tini:tend,:,yini:yend,:].transpose()
+                    psfct = psfc[tini:tend,yini:yend,:].transpose()
+                    hgtt = hgt[yini:yend,:].transpose()
+                    tempt = temp[tini:tend,:,yini:yend,:].transpose()
+                    qvt = qv[tini:tend,:,yini:yend,:].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()
+
+# Finishing split along y-axis
+            if dimyfrac != dimy and yend != dimy-2:
+                yini = yend + 0
+                yend = dimy
+                if not ysaid:
+                    print '    finishing yaxis-splitting:', yini, ', ', yend
+                    ysaid = True
+
+                Ndy = len(range(yini,yend))
+                if gen.searchInlist(notCHK, vn):
+                    if vn == 'WRFght': 
+                        varin = geop[tini:tend,:,yini:yend,:]
+                        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[tini:tend,:,yini:yend,:]/(1.+qv[tini:tend,:,yini:yend,:])
+                        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[tini:tend,:,yini:yend,:]
+                        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']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering x-wind variable !!'
+                        print '    from:', ovarin.shape, 'to', (dimt, dimz, dimy, dimx)
+                        if isFR64:
+                            varin = np.zeros((Ndt, dimz, Ndy, dimx), dtype=np.float64)
+                        else:
+                            varin = np.zeros((Ndt, dimz, Ndy, dimx), dtype=np.float)
+                        varin[tini:tend,:,yini:yend,:] =                             \
+                          0.5*(ovarin[tini:tend,:,yini:yend,0:dimx] +                \
+                          ovarin[tini:tend,:,yini:yend,1:dimx+1])
+
+                        # Not pro, but less memory problems!
+                        #for it in range(dimt):
+                        #    varin[it,:,::] = 0.5*(ovarin[it,:,:,0:dimx] + ovarin[it,:,:,1:dimx+1])
+                        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']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering x-wind variable & ' +  \
+                          'Earth rotating!!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        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,:,yini:yend,:] =                               \
+                              0.5*(ovarin[it,:,yini:yend,0:dimx] +                   \
+                              ovarin[it,:,yini:yend,1:dimx+1])
+                            varin02[it,:,yini:yend,:] =                              \
+                              0.5*(ovarin2[it,:,yini:yend,:] +                       \
+                              ovarin2[it,:,yini+1:yend+1,:])
+                        for iz in range(dimz):
+                            varin[tini:tend,iz,yini:yend,:] =                        \
+                              varin0[tini:tend,iz,yini:yend,:]*                      \
+                              ocosa[tini:tend,yini:yend,:] -                         \
+                              varin02[tini:tend,iz,yini:yend,:]*                     \
+                              osina[tini:tend,yini:yend,:]
+                        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 == 'WRFv':
+                        ovarin = onc.variables['V']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering y-wind variable !!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        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,:,yini:yend,:] = 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']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering y-wind variable & ' +  \
+                          'Earth rotating!!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        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,:,yini:yend,:] =                               \
+                              0.5*(ovarin[it,:,yini:yend,0:dimx] +                   \
+                              ovarin[it,:,yini:yend,1:dimx+1])
+                            varin02[it,:,yini:yend,:] =                              \
+                              0.5*(ovarin2[it,:,yini:yend,:] +                       \
+                              ovarin2[it,:,yini+1:yend+1,:])
+                        for iz in range(dimz):
+                            varin[tini:tend,iz,yini:yend,:] =                        \
+                              varin0[tini:tend,iz,yini:yend,:]*                      \
+                              osina[tini:tend,yini:yend,:] +                         \
+                              varin02[tini:tend,iz,yini:yend,:]*                     \
+                              ocosa[tini:tend,yini:yend,:]
+                        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 == 'WRFw':
+                        ovarin = onc.variables['W']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering z-wind variable !!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        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,:,yini:yend,:] =                                \
+                              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[tini:tend,:,yini:yend,:]
+                        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[tini:tend,:,yini:yend,:]
+                    # 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:
+                    varint = varin.transpose()
+                    prest = pres[tini:tend,:,yini:yend,:].transpose()
+                    psfct = psfc[tini:tend,yini:yend,:].transpose()
+                    hgtt = hgt[yini:yend,:].transpose()
+                    tempt = temp[tini:tend,:,yini:yend,:].transpose()
+                    qvt = qv[tini:tend,:,yini:yend,:].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()
+
+        # Finishing dimtfrac
+        if dimtfrac != dimt and tend != dimt-2:
+            if not tsaid:
+                print '    finishing time-splitting:', tend, ', ', dimt
+                tsaid = True
+            tini = tend
+            tend = dimt
+            Ndt = len(range(tini,tend))
+            for iyy in range(0,dimyfrac*int(dimy/dimyfrac),dimyfrac):
+                yini = iyy
+                yend = iyy + dimyfrac
+                Ndy = len(range(yini,yend))
+
+                if gen.searchInlist(notCHK, vn):
+                    if vn == 'WRFght': 
+                        varin = geop[tini:tend,:,yini:yend,:]
+                        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[tini:tend,:,yini:yend,:]/(1.+qv[tini:tend,:,yini:yend,:])
+                        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[tini:tend,:,yini:yend,:]
+                        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']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering x-wind variable !!'
+                        print '    from:', ovarin.shape, 'to', (dimt, dimz, dimy, dimx)
+                        if isFR64:
+                            varin = np.zeros((Ndt, dimz, Ndy, dimx), dtype=np.float64)
+                        else:
+                            varin = np.zeros((Ndt, dimz, Ndy, dimx), dtype=np.float)
+                        varin[tini:tend,:,yini:yend,:] =                             \
+                          0.5*(ovarin[tini:tend,:,yini:yend,0:dimx] +                \
+                          ovarin[tini:tend,:,yini:yend,1:dimx+1])
+
+                        # Not pro, but less memory problems!
+                        #for it in range(dimt):
+                        #    varin[it,:,::] = 0.5*(ovarin[it,:,:,0:dimx] + ovarin[it,:,:,1:dimx+1])
+                        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']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering x-wind variable & ' +  \
+                          'Earth rotating!!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        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,:,yini:yend,:] =                               \
+                              0.5*(ovarin[it,:,yini:yend,0:dimx] +                   \
+                              ovarin[it,:,yini:yend,1:dimx+1])
+                            varin02[it,:,yini:yend,:] =                              \
+                              0.5*(ovarin2[it,:,yini:yend,:] +                       \
+                              ovarin2[it,:,yini+1:yend+1,:])
+                        for iz in range(dimz):
+                            varin[tini:tend,iz,yini:yend,:] =                        \
+                              varin0[tini:tend,iz,yini:yend,:]*                      \
+                              ocosa[tini:tend,yini:yend,:] -                         \
+                              varin02[tini:tend,iz,yini:yend,:]*                     \
+                              osina[tini:tend,yini:yend,:]
+                        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 == 'WRFv':
+                        ovarin = onc.variables['V']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering y-wind variable !!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        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,:,yini:yend,:] = 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']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering y-wind variable & ' +  \
+                          'Earth rotating!!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        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,:,yini:yend,:] =                               \
+                              0.5*(ovarin[it,:,yini:yend,0:dimx] +                   \
+                              ovarin[it,:,yini:yend,1:dimx+1])
+                            varin02[it,:,yini:yend,:] =                              \
+                              0.5*(ovarin2[it,:,yini:yend,:] +                       \
+                              ovarin2[it,:,yini+1:yend+1,:])
+                        for iz in range(dimz):
+                            varin[tini:tend,iz,yini:yend,:] =                        \
+                              varin0[tini:tend,iz,yini:yend,:]*                      \
+                              osina[tini:tend,yini:yend,:] +                         \
+                              varin02[tini:tend,iz,yini:yend,:]*                     \
+                              ocosa[tini:tend,yini:yend,:]
+                        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 == 'WRFw':
+                        ovarin = onc.variables['W']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering z-wind variable !!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        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,:,yini:yend,:] =                                \
+                              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[tini:tend,:,yini:yend,:]
+                        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[tini:tend,:,yini:yend,:]
+                    # 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:
+                    varint = varin.transpose()
+                    prest = pres[tini:tend,:,yini:yend,:].transpose()
+                    psfct = psfc[tini:tend,yini:yend,:].transpose()
+                    hgtt = hgt[yini:yend,:].transpose()
+                    tempt = temp[tini:tend,:,yini:yend,:].transpose()
+                    qvt = qv[tini:tend,:,yini:yend,:].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()
+
+# Finishing split along y-axis
+            if dimyfrac != dimy and yend != dimy-2:
+                yini = yend + 0
+                yend = dimy
+                if not ysaid:
+                    print '    finishing yaxis-splitting:', yini, ', ', yend
+                    ysaid = True
+
+                Ndy = len(range(yini,yend))
+                if gen.searchInlist(notCHK, vn):
+                    if vn == 'WRFght': 
+                        varin = geop[tini:tend,:,yini:yend,:]
+                        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[tini:tend,:,yini:yend,:]/(1.+qv[tini:tend,:,yini:yend,:])
+                        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[tini:tend,:,yini:yend,:]
+                        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']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering x-wind variable !!'
+                        print '    from:', ovarin.shape, 'to', (dimt, dimz, dimy, dimx)
+                        if isFR64:
+                            varin = np.zeros((Ndt, dimz, Ndy, dimx), dtype=np.float64)
+                        else:
+                            varin = np.zeros((Ndt, dimz, Ndy, dimx), dtype=np.float)
+                        varin[tini:tend,:,yini:yend,:] =                             \
+                          0.5*(ovarin[tini:tend,:,yini:yend,0:dimx] +                \
+                          ovarin[tini:tend,:,yini:yend,1:dimx+1])
+
+                        # Not pro, but less memory problems!
+                        #for it in range(dimt):
+                        #    varin[it,:,::] = 0.5*(ovarin[it,:,:,0:dimx] + ovarin[it,:,:,1:dimx+1])
+                        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']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering x-wind variable & ' +  \
+                          'Earth rotating!!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        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,:,yini:yend,:] =                               \
+                              0.5*(ovarin[it,:,yini:yend,0:dimx] +                   \
+                              ovarin[it,:,yini:yend,1:dimx+1])
+                            varin02[it,:,yini:yend,:] =                              \
+                              0.5*(ovarin2[it,:,yini:yend,:] +                       \
+                              ovarin2[it,:,yini+1:yend+1,:])
+                        for iz in range(dimz):
+                            varin[tini:tend,iz,yini:yend,:] =                        \
+                              varin0[tini:tend,iz,yini:yend,:]*                      \
+                              ocosa[tini:tend,yini:yend,:] -                         \
+                              varin02[tini:tend,iz,yini:yend,:]*                     \
+                              osina[tini:tend,yini:yend,:]
+                        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 == 'WRFv':
+                        ovarin = onc.variables['V']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering y-wind variable !!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        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,:,yini:yend,:] = 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']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering y-wind variable & ' +  \
+                          'Earth rotating!!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        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,:,yini:yend,:] =                               \
+                              0.5*(ovarin[it,:,yini:yend,0:dimx] +                   \
+                              ovarin[it,:,yini:yend,1:dimx+1])
+                            varin02[it,:,yini:yend,:] =                              \
+                              0.5*(ovarin2[it,:,yini:yend,:] +                       \
+                              ovarin2[it,:,yini+1:yend+1,:])
+                        for iz in range(dimz):
+                            varin[tini:tend,iz,yini:yend,:] =                        \
+                              varin0[tini:tend,iz,yini:yend,:]*                      \
+                              osina[tini:tend,yini:yend,:] +                         \
+                              varin02[tini:tend,iz,yini:yend,:]*                     \
+                              ocosa[tini:tend,yini:yend,:]
+                        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 == 'WRFw':
+                        ovarin = onc.variables['W']
+                        print infmsg
+                        print '  ' + fname + ': De-staggering z-wind variable !!'
+                        print '    from:', ovarin.shape, 'to', (dimt,dimz,dimy,dimx)
+                        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,:,yini:yend,:] =                                \
+                              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[tini:tend,:,yini:yend,:]
+                        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[tini:tend,:,yini:yend,:]
+                    # 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:
+                    varint = varin.transpose()
+                    prest = pres[tini:tend,:,yini:yend,:].transpose()
+                    psfct = psfc[tini:tend,yini:yend,:].transpose()
+                    hgtt = hgt[yini:yend,:].transpose()
+                    tempt = temp[tini:tend,:,yini:yend,:].transpose()
+                    qvt = qv[tini:tend,:,yini:yend,:].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
+
+    onewnc.setncattr('author', 'L. Fita')
+    newattr = set_attributek(onewnc, 'institution', unicode('Laboratoire de M' +     \
+      unichr(233) + 't' + unichr(233) + 'orologie Dynamique'), 'U')
+    onewnc.setncattr('university', 'Pierre Marie Curie - Jussieu')
+    onewnc.setncattr('center', 'Centre National de Recherches Scientifiques')
+    onewnc.setncattr('city', 'Paris')
+    onewnc.setncattr('country', 'France')
+    onewnc.setncattr('script', 'nc_var_tools.py')
+    onewnc.setncattr('function', 'pinterp')
+    onewnc.setncattr('version', '1.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_simple(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]
         [interplevs]: ':' separated list of pressure values to interpolate [Pa] (decreassing in Pa)
@@ -18774 +20085 @@
     return
 
+
 def get_Variables(values, ncfile, variables):
     """ Function to get a list of variables from an existing file