# Python script to comput diagnostics # From L. Fita work in different places: CCRC (Australia), LMD (France) # More information at: http://www.xn--llusfb-5va.cat/python/PyNCplot # # pyNCplot and its component nc_var.py comes with ABSOLUTELY NO WARRANTY. # This work is licendes under a Creative Commons # Attribution-ShareAlike 4.0 International License (http://creativecommons.org/licenses/by-sa/4.0) # # L. Fita, CIMA. CONICET-UBA, CNRS UMI-IFAECI, C.A. Buenos Aires, Argentina # File diagnostics.inf provides the combination of variables to get the desired diagnostic # To be used with module_ForDiagnostics.F90, module_ForDiagnosticsVars.F90, module_generic.F90 # foudre: f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90 >& run_f2py.log # ciclad: f2py --f90flags="-fPIC" --f90exec=/usr/bin/gfortran -L/opt/canopy-1.3.0/Canopy_64bit/System/lib/ -L/usr/lib64/ -L/opt/canopy-1.3.0/Canopy_64bit/System/lib/ -m module_ForDiagnostics -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90 >& run_f2py.log ## e.g. # diagnostics.py -d 'Time@WRFtime,bottom_top@ZNU,south_north@XLAT,west_east@XLONG' -v 'clt|CLDFRA,cllmh|CLDFRA@WRFp,RAINTOT|RAINC@RAINNC@RAINSH@XTIME' -f WRF_LMDZ/NPv31/wrfout_d01_1980-03-01_00:00:00 ## e.g. # diagnostics.py -f /home/lluis/PY/diagnostics.inf -d variable_combo -v WRFprc # Available general pupose diagnostics (model independent) providing (varv1, varv2, ..., dimns, dimvns) # compute_accum: Function to compute the accumulation of a variable # compute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following # newmicro.F90 from LMDZ compute_clt(cldfra, pres, dimns, dimvns) # compute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ # compute_clivi: Function to compute cloud-ice water path (clivi) # compute_clwvl: Function to compute condensed water path (clwvl) # compute_deaccum: Function to compute the deaccumulation of a variable # compute_mslp: Function to compute mslp: mean sea level pressure following p_interp.F90 from WRF # compute_OMEGAw: Function to transform OMEGA [Pas-1] to velocities [ms-1] # compute_prw: Function to compute water vapour path (prw) # compute_range_faces: Function to compute faces [uphill, valley, downhill] of sections of a mountain # range, along a given face # compute_rh: Function to compute relative humidity following 'Tetens' equation (T,P) ...' # compute_td: Function to compute the dew point temperature # compute_turbulence: Function to compute the rubulence term of the Taylor's decomposition ...' # compute_wds: Function to compute the wind direction # compute_wss: Function to compute the wind speed # compute_WRFuava: Function to compute geographical rotated WRF 3D winds # compute_WRFuasvas: Fucntion to compute geographical rotated WRF 2-meter winds # derivate_centered: Function to compute the centered derivate of a given field # def Forcompute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ via Fortran subroutine # Forcompute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ via a Fortran module # Forcompute_psl_ptarget: Function to compute the sea-level pressure following target_pressure value found in `p_interp.F' # Others just providing variable values # var_cllmh: Fcuntion to compute cllmh on a 1D column # var_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ using 1D vertical column values # var_mslp: Fcuntion to compute mean sea-level pressure # var_virtualTemp: This function returns virtual temperature in K, # var_WRFtime: Function to copmute CFtimes from WRFtime variable # rotational_z: z-component of the rotatinoal of horizontal vectorial field # turbulence_var: Function to compute the Taylor's decomposition turbulence term from a a given variable # timeoverthres: When a given variable [varname] overpass a given [value]. Being [CFvarn] the name of the diagnostics in # variables_values.dat # timemax ([varname], time). When a given variable [varname] got its maximum from optparse import OptionParser import numpy as np import numpy.ma as ma from netCDF4 import Dataset as NetCDFFile import os import re import nc_var_tools as ncvar import generic_tools as gen import datetime as dtime import module_ForDiag as fdin import module_ForDef as fdef import diag_tools as diag main = 'diagnostics.py' errormsg = 'ERROR -- error -- ERROR -- error' warnmsg = 'WARNING -- warning -- WARNING -- warning' # Constants grav = 9.81 ####### ###### ##### #### ### ## # comboinf="\nIF -d 'variable_combo', provides information of the combination to obtain -v [varn] with the ASCII file with the combinations as -f [combofile]" parser = OptionParser() parser.add_option("-f", "--netCDF_file", dest="ncfile", help="file to use", metavar="FILE") parser.add_option("-d", "--dimensions", dest="dimns", help="[dimtn]@[dtvn],[dimzn]@[dzvn],[...,[dimxn]@[dxvn]], ',' list with the couples [dimDn]@[dDvn], [dimDn], name of the dimension D and name of the variable [dDvn] with the values of the dimension ('WRFtime', for WRF time copmutation). NOTE: same order as in file!!!!" + comboinf, metavar="LABELS") parser.add_option("-v", "--variables", dest="varns", help=" [varn1]|[var11]@[...[varN1]],[...,[varnM]|[var1M]@[...[varLM]]] ',' list of variables to compute [varnK] and its necessary ones [var1K]...[varPK]", metavar="VALUES") (opts, args) = parser.parse_args() ####### ####### ## MAIN ####### availdiags = ['ACRAINTOT', 'accum', 'clt', 'cllmh', 'convini', 'deaccum', 'fog_K84', \ 'fog_RUC', 'front_R04', 'frontogenesis', 'gradient2Dh', 'rhs_tas_tds', 'LMDZrh', \ 'mslp', 'OMEGAw', 'RAINTOT', 'range_faces', \ 'rvors', 'td', 'timemax', 'timeoverthres', 'turbulence', 'tws', 'uavaFROMwswd', \ 'WRFbnds', 'WRFcape_afwa', 'WRFclivi', 'WRFclwvi', 'WRF_denszint', 'WRFgeop', \ 'WRFmrso', 'WRFmrsos', 'WRFpotevap_orPM', 'WRFp', 'WRFpsl_ecmwf', \ 'WRFpsl_ptarget', 'WRFrvors', 'WRFslw', 'ws', 'wds', 'wss', 'WRFheight', \ 'WRFheightrel', 'WRFtda', 'WRFtdas', 'WRFtws', 'WRFua', 'WRFva', 'WRFzwind', \ 'WRFzwind_log', 'WRFzwindMO', 'zmlagen', 'zmlagenUWsnd'] methods = ['accum', 'deaccum'] # Variables not to check NONcheckingvars = ['accum', 'cllmh', 'deaccum', 'face', 'LONLATdxdy', 'LONLATdx', \ 'LONLATdy', 'params', 'reglonlatbnds', 'TSrhs', 'TStd', 'TSwds', 'TSwss', \ 'UNua', 'UNva', 'UNwa', \ 'WRFbils', 'WRFbnds', \ 'WRFclivi', 'WRFclwvi', 'WRFdens', 'WRFdx', 'WRFdxdy', 'WRFdxdywps', 'WRFdy', \ 'WRFdz', 'WRFgeop', 'WRFp', 'WRFtd', \ 'WRFpos', 'WRFprc', 'WRFprls', 'WRFrh', 'LMDZrh', 'LMDZrhs', \ 'WRFrhs', 'WRFrvors', \ 'WRFt', 'WRFtime', 'WRFua', 'WRFva', 'WRFwds', 'WRFwss', 'WRFheight', 'WRFz', \ 'WRFzg'] # diagnostics not to check their dependeny NONcheckdepvars = ['accum', 'deaccum', 'timeoverthres', 'WRF_denszint', \ 'WRFzwind_log', 'WRFzwind', 'WRFzwindMO'] NONchkvardims = ['WRFtime'] ofile = 'diagnostics.nc' dimns = opts.dimns varns = opts.varns # Special method. knowing variable combination ## if opts.dimns == 'variable_combo': print warnmsg print ' ' + main + ': knowing variable combination !!!' combination = variable_combo(opts.varns,opts.ncfile) print ' COMBO: ' + combination quit(-1) if opts.ncfile is None: print errormsg print ' ' + main + ": No file provided !!" print ' is mandatory to provide a file -f [filename]' quit(-1) if opts.dimns is None: print errormsg print ' ' + main + ": No description of dimensions are provided !!" print ' is mandatory to provide description of dimensions as -d [dimn]@[vardimname],... ' quit(-1) if opts.varns is None: print errormsg print ' ' + main + ": No variable to diagnose is provided !!" print ' is mandatory to provide a variable to diagnose as -v [diagn]|[varn1]@... ' quit(-1) if not os.path.isfile(opts.ncfile): print errormsg print ' ' + main + ": file '" + opts.ncfile + "' does not exist !!" quit(-1) ncobj = NetCDFFile(opts.ncfile, 'r') # Looking for specific variables that might be use in more than one diagnostic WRFgeop_compute = False WRFp_compute = False WRFt_compute = False WRFrh_compute = False WRFght_compute = False WRFdens_compute = False WRFpos_compute = False WRFtime_compute = False WRFz_compute = False WRFdx_compute = False WRFdy_compute = False WRFdz_compute = False WRFdxdy_compute = False WRFdxdywps_compute = False LONLATdxdy_compute = False LONLATdx_compute = False LONLATdy_compute = False UNua_compute = False UNva_compute = False UNwa_compute = False # File creation newnc = NetCDFFile(ofile,'w') # dimensions dimvalues = dimns.split(',') dnames = [] dvnames = [] for dimval in dimvalues: dn = dimval.split('@')[0] dnv = dimval.split('@')[1] dnames.append(dn) dvnames.append(dnv) # Is there any dimension-variable which should be computed? if dnv == 'WRFgeop':WRFgeop_compute = True if dnv == 'WRFp': WRFp_compute = True if dnv == 'WRFt': WRFt_compute = True if dnv == 'WRFrh': WRFrh_compute = True if dnv == 'WRFght': WRFght_compute = True if dnv == 'WRFdens': WRFdens_compute = True if dnv == 'WRFpos': WRFpos_compute = True if dnv == 'WRFtime': WRFtime_compute = True if dnv == 'WRFz':WRFz_compute = True if dnv == 'WRFdx':WRFdx_compute = True if dnv == 'WRFdy':WRFdy_compute = True if dnv == 'WRFdz':WRFdz_compute = True if dnv == 'WRFdxdy':WRFdxdy_compute = True if dnv == 'WRFdxdywps':WRFdxdywps_compute = True if dnv == 'LONLATdxdy':LONLATdxdy_compute = True if dnv == 'LONLATdx':LONLATdx_compute = True if dnv == 'LONLATdy':LONLATdy_compute = True if dnv[0:4] == 'UNua': UNua_compute = True vUnua = dnv.split(',')[1] if dnv[0:4] == 'UNva': UNva_compute = True vUnva = dnv.split(',')[1] if dnv[0:4] == 'UNwa': UNwa_compute = True vUnwa = dnv.split(',')[1] # diagnostics to compute diags = varns.split(',') Ndiags = len(diags) for idiag in range(Ndiags): if diags[idiag].split('|')[1].find('@') == -1: depvars = diags[idiag].split('|')[1] if depvars == 'WRFgeop':WRFgeop_compute = True if depvars == 'WRFp': WRFp_compute = True if depvars == 'WRFt': WRFt_compute = True if depvars == 'WRFrh': WRFrh_compute = True if depvars == 'WRFght': WRFght_compute = True if depvars == 'WRFdens': WRFdens_compute = True if depvars == 'WRFpos': WRFpos_compute = True if depvars == 'WRFtime': WRFtime_compute = True if depvars == 'WRFz': WRFz_compute = True if depvars == 'WRFdx': WRFdx_compute = True if depvars == 'WRFdy': WRFdy_compute = True if depvars == 'WRFdz': WRFdz_compute = True if depvars == 'WRFdxdy': WRFdxdy_compute = True if depvars == 'WRFdxdywps': WRFdxdywps_compute = True if depvars == 'LONLATdxdy': LONLATdxdy_compute = True if depvars == 'LONLATdx': LONLATdx_compute = True if depvars == 'LONLATdy': LONLATdy_compute = True if depvars[0:4] == 'UNua': UNua_compute = True vUnua = dnv.split(',')[1] if depvars[0:4] == 'UNva': UNva_compute = True vUnva = dnv.split(',')[1] if depvars[0:4] == 'UNwa': UNwa_compute = True vUnwa = dnv.split(',')[1] else: depvars = diags[idiag].split('|')[1].split('@') if gen.searchInlist(depvars, 'WRFgeop'): WRFgeop_compute = True if gen.searchInlist(depvars, 'WRFp'): WRFp_compute = True if gen.searchInlist(depvars, 'WRFt'): WRFt_compute = True if gen.searchInlist(depvars, 'WRFrh'): WRFrh_compute = True if gen.searchInlist(depvars, 'WRFght'): WRFght_compute = True if gen.searchInlist(depvars, 'WRFdens'): WRFdens_compute = True if gen.searchInlist(depvars, 'WRFpos'): WRFpos_compute = True if gen.searchInlist(depvars, 'WRFtime'): WRFtime_compute = True if gen.searchInlist(depvars, 'WRFz'): WRFz_compute = True if gen.searchInlist(depvars, 'WRFdx'): WRFdx_compute = True if gen.searchInlist(depvars, 'WRFdy'): WRFdy_compute = True if gen.searchInlist(depvars, 'WRFdz'): WRFdz_compute = True if gen.searchInlist(depvars, 'WRFdxdy'): WRFdxdy_compute = True if gen.searchInlist(depvars, 'WRFdxdywps'): WRFdxdywps_compute = True if gen.searchInlist(depvars, 'LONLATdxdy'): LONLATdxdy_compute = True if gen.searchInlist(depvars, 'LONLATdx'): LONLATdx_compute = True if gen.searchInlist(depvars, 'LONLATdy'): LONLATdy_compute = True if gen.searchInlist_Strsec(depvars, 0, 3, 'UNua'): UNua_compute = True vals, ind = gen.search_sec_list(depvars,'UNua') dnv = depvars[ind[0]] vUNua = dnv.split(':')[1] if gen.searchInlist_Strsec(depvars, 0, 3, 'UNva'): UNva_compute = True vals, ind = gen.search_sec_list(depvars,'UNva') dnv = depvars[ind[0]] vUNva = dnv.split(':')[1] if gen.searchInlist_Strsec(depvars, 0, 3, 'UNwa'): UNwa_compute = True vals, ind = gen.search_sec_list(depvars,'UNwa') dnv = depvars[ind[0]] vUNwa = dnv.split(':')[1] # Dictionary with the new computed variables to be able to add them dictcompvars = {} if WRFgeop_compute: print ' ' + main + ': Retrieving geopotential value from WRF as PH + PHB' dimv = ncobj.variables['PH'].shape WRFgeop = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:] # Attributes of the variable Vvals = gen.variables_values('WRFgeop') dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} if WRFp_compute: print ' ' + main + ': Retrieving pressure value from WRF as P + PB' dimv = ncobj.variables['P'].shape WRFp = ncobj.variables['P'][:] + ncobj.variables['PB'][:] # Attributes of the variable Vvals = gen.variables_values('WRFp') dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} if WRFght_compute: print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...' WRFght = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:] # Attributes of the variable Vvals = gen.variables_values('WRFght') dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} if WRFrh_compute: print ' ' + main + ": computing relative humidity from WRF as 'Tetens'" + \ ' equation (T,P) ...' p0=100000. p=ncobj.variables['P'][:] + ncobj.variables['PB'][:] tk = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.) qv = ncobj.variables['QVAPOR'][:] data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65)) data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1) WRFrh = qv/data2 # Attributes of the variable Vvals = gen.variables_values('WRFrh') dictcompvars['WRFrh'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} if WRFt_compute: print ' ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...' p0=100000. p=ncobj.variables['P'][:] + ncobj.variables['PB'][:] WRFt = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.) # Attributes of the variable Vvals = gen.variables_values('WRFt') dictcompvars['WRFt'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} if WRFdens_compute: print ' ' + main + ': computing air density from WRF as ((MU + MUB) * ' + \ 'DNW)/g ...' # Just we need in in absolute values: Size of the central grid cell ## dxval = ncobj.getncattr('DX') ## dyval = ncobj.getncattr('DY') ## mapfac = ncobj.variables['MAPFAC_M'][:] ## area = dxval*dyval*mapfac mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:]) dnw = ncobj.variables['DNW'][:] WRFdens = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]), \ dtype=np.float) levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float) for it in range(mu.shape[0]): for iz in range(dnw.shape[1]): levval.fill(np.abs(dnw[it,iz])) WRFdens[it,iz,:,:] = levval WRFdens[it,iz,:,:] = mu[it,:,:]*WRFdens[it,iz,:,:]/grav # Attributes of the variable Vvals = gen.variables_values('WRFdens') dictcompvars['WRFdens'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} if WRFpos_compute: # WRF positions from the lowest-leftest corner of the matrix print ' ' + main + ': computing position from MAPFAC_M as sqrt(DY*j**2 + ' + \ 'DX*x**2)*MAPFAC_M ...' mapfac = ncobj.variables['MAPFAC_M'][:] distx = np.float(ncobj.getncattr('DX')) disty = np.float(ncobj.getncattr('DY')) print 'distx:',distx,'disty:',disty dx = mapfac.shape[2] dy = mapfac.shape[1] dt = mapfac.shape[0] WRFpos = np.zeros((dt, dy, dx), dtype=np.float) for i in range(1,dx): WRFpos[0,0,i] = distx*i/mapfac[0,0,i] for j in range(1,dy): i=0 WRFpos[0,j,i] = WRFpos[0,j-1,i] + disty/mapfac[0,j,i] for i in range(1,dx): # WRFpos[0,j,i] = np.sqrt((disty*j)**2. + (distx*i)**2.)/mapfac[0,j,i] # WRFpos[0,j,i] = np.sqrt((disty*j)**2. + (distx*i)**2.) WRFpos[0,j,i] = WRFpos[0,j,i-1] + distx/mapfac[0,j,i] for it in range(1,dt): WRFpos[it,:,:] = WRFpos[0,:,:] if WRFtime_compute: print ' ' + main + ': computing time from WRF as CFtime(Times) ...' refdate='19491201000000' tunitsval='minutes' timeobj = ncobj.variables['Times'] timewrfv = timeobj[:] yrref=refdate[0:4] monref=refdate[4:6] dayref=refdate[6:8] horref=refdate[8:10] minref=refdate[10:12] secref=refdate[12:14] refdateS = yrref + '-' + monref + '-' + dayref + ' ' + horref + ':' + minref + \ ':' + secref if len(timeobj.shape) == 2: dt = timeobj.shape[0] else: dt = 1 WRFtime = np.zeros((dt), dtype=np.float) if len(timeobj.shape) == 2: for it in range(dt): wrfdates = gen.datetimeStr_conversion(timewrfv[it,:],'WRFdatetime', \ 'matYmdHMS') WRFtime[it] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval) else: wrfdates = gen.datetimeStr_conversion(timewrfv[:],'WRFdatetime', \ 'matYmdHMS') WRFtime[0] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval) tunits = tunitsval + ' since ' + refdateS # Attributes of the variable dictcompvars['WRFtime'] = {'name': 'time', 'standard_name': 'time', \ 'long_name': 'time', 'units': tunits, 'calendar': 'gregorian'} if WRFz_compute: print ' ' + main + ': Retrieving z: height above surface value from WRF as ' + \ 'unstagger(PH + PHB)/9.8-hgt' dimv = ncobj.variables['PH'].shape WRFzg = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/9.8 unzgd = (dimv[0], dimv[1]-1, dimv[2], dimv[3]) unzg = np.zeros(unzgd, dtype=np.float) unzg = 0.5*(WRFzg[:,0:dimv[1]-1,:,:] + WRFzg[:,1:dimv[1],:,:]) WRFz = np.zeros(unzgd, dtype=np.float) for iz in range(dimv[1]-1): WRFz[:,iz,:,:] = unzg[:,iz,:,:] - ncobj.variables['HGT'][:] # Attributes of the variable Vvals = gen.variables_values('WRFz') dictcompvars['WRFz'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} if WRFdxdy_compute or WRFdx_compute or WRFdy_compute: print ' ' + main + ': Retrieving dxdy: real distance between grid points ' + \ 'from WRF as dx=(XLONG(i+1)-XLONG(i))*DX/MAPFAC_M, dy=(XLAT(j+1)-XLAT(i))*DY/'+\ 'MAPFAC_M, ds=sqrt(dx**2+dy**2)' dimv = ncobj.variables['XLONG'].shape WRFlon = ncobj.variables['XLONG'][0,:,:] WRFlat = ncobj.variables['XLAT'][0,:,:] WRFmapfac_m = ncobj.variables['MAPFAC_M'][0,:,:] DX = ncobj.DX DY = ncobj.DY dimx = dimv[2] dimy = dimv[1] WRFdx = np.zeros((dimy,dimx), dtype=np.float) WRFdy = np.zeros((dimy,dimx), dtype=np.float) WRFdx[:,0:dimx-1]=(WRFlon[:,1:dimx]-WRFlon[:,0:dimx-1])*DX/WRFmapfac_m[:,0:dimx-1] WRFdy[0:dimy-1,:]=(WRFlat[1:dimy,:]-WRFlat[0:dimy-1,:])*DY/WRFmapfac_m[0:dimy-1,:] WRFds = np.sqrt(WRFdx**2 + WRFdy**2) # Attributes of the variable Vvals = gen.variables_values('WRFdx') dictcompvars['WRFdx'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} Vvals = gen.variables_values('WRFdy') dictcompvars['WRFdy'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} Vvals = gen.variables_values('WRFds') dictcompvars['WRFds'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} if WRFdxdywps_compute: print ' ' + main + ': Retrieving dxdy: real distance between grid points ' + \ 'from wpsWRF as dx=(XLONG_M(i+1)-XLONG_M(i))*DX/MAPFAC_M, ' + \ 'dy=(XLAT_M(j+1)-XLAT_M(i))*DY/MAPFAC_M, ds=sqrt(dx**2+dy**2)' dimv = ncobj.variables['XLONG_M'].shape WRFlon = ncobj.variables['XLONG_M'][0,:,:] WRFlat = ncobj.variables['XLAT_M'][0,:,:] WRFmapfac_m = ncobj.variables['MAPFAC_M'][0,:,:] DX = ncobj.DX DY = ncobj.DY dimx = dimv[2] dimy = dimv[1] WRFdx = np.zeros((dimy,dimx), dtype=np.float) WRFdy = np.zeros((dimy,dimx), dtype=np.float) WRFdx[:,0:dimx-1]=(WRFlon[:,1:dimx]-WRFlon[:,0:dimx-1])*DX/WRFmapfac_m[:,0:dimx-1] WRFdy[0:dimy-1,:]=(WRFlat[1:dimy,:]-WRFlat[0:dimy-1,:])*DY/WRFmapfac_m[0:dimy-1,:] WRFds = np.sqrt(WRFdx**2 + WRFdy**2) # Attributes of the variable Vvals = gen.variables_values('WRFdx') dictcompvars['WRFdx'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} Vvals = gen.variables_values('WRFdy') dictcompvars['WRFdy'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} Vvals = gen.variables_values('WRFds') dictcompvars['WRFds'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} if WRFdz_compute: print ' ' + main + ': Retrieving dz: real distance between grid points ' + \ 'from WRF as dz=PHB(k+1)+PH(k+1)-(PHB(k)+PH(k))' PH = ncobj.variables['PH'][:] PHB = ncobj.variables['PHB'][:] dimv = ncobj.variables['PH'].shape dimx = dimv[3] dimy = dimv[2] dimz = dimv[1] dimt = dimv[0] WRFdz = np.zeros((dimt,dimz,dimy,dimx), dtype=np.float) WRFdz[:,0:dimz-1,:,:]=PH[:,1:dimz,:,:]+PHB[:,1:dimz,:,:]-(PH[:,0:dimz-1,:,:]+ \ PHB[:,0:dimz-1,:,:]) # Attributes of the variable Vvals = gen.variables_values('WRFdz') dictcompvars['WRFdz'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} if LONLATdxdy_compute or LONLATdx_compute or LONLATdy_compute : print ' ' + main + ': Retrieving dxdy: real distance between grid points ' + \ 'from a regular lonlat projection as dx=(lon[i+1]-lon[i])*raddeg*Rearth*' + \ 'cos(abs(lat[i])); dy=(lat[j+1]-lat[i])*raddeg*Rearth; ds=sqrt(dx**2+dy**2); '+\ 'raddeg = pi/180; Rearth=6370.0e03' dimv = ncobj.variables['lon'].shape lon = ncobj.variables['lon'][:] lat = ncobj.variables['lat'][:] WRFlon, WRFlat = gen.lonlat2D(lon,lat) dimx = WRFlon.shape[1] dimy = WRFlon.shape[0] WRFdx = np.zeros((dimy,dimx), dtype=np.float) WRFdy = np.zeros((dimy,dimx), dtype=np.float) raddeg = np.pi/180. Rearth = fdef.module_definitions.earthradii WRFdx[:,0:dimx-1]=(WRFlon[:,1:dimx]-WRFlon[:,0:dimx-1])*raddeg*Rearth* \ np.cos(np.abs(WRFlat[:,0:dimx-1]*raddeg)) WRFdy[0:dimy-1,:]=(WRFlat[1:dimy,:]-WRFlat[0:dimy-1,:])*raddeg*Rearth WRFds = np.sqrt(WRFdx**2 + WRFdy**2) # Attributes of the variable Vvals = gen.variables_values('WRFdx') dictcompvars['WRFdx'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} Vvals = gen.variables_values('WRFdy') dictcompvars['WRFdy'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} Vvals = gen.variables_values('WRFds') dictcompvars['WRFds'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} if UNua_compute: print ' ' + main + ": un-staggering '" + vUNua + "': as 0.5*(ua[0:dimx-1]+" + \ "ua[1:dimx])" vunua = ncobj.variables[vUNua][:] dimv = ncobj.variables[vUNua].shape dimx = dimv[3] dimy = dimv[2] dimz = dimv[1] dimt = dimv[0] undimx = 'unx' unua = np.zeros((dimt,dimz,dimy,dimx-1), dtype=np.float) unua[...,0:dimx-1] = 0.5*(vunua[...,1:dimx]+vunua[...,0:dimx-1]) # Attributes of the variable Vvals = gen.variables_values('ua') dictcompvars['unua'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} if UNva_compute: print ' ' + main + ": un-staggering '" + vUNva + "': as 0.5*(va[0:dimy-1]+" + \ "va[1:dimy])" vunva = ncobj.variables[vUNva][:] dimv = ncobj.variables[vUNva].shape dimx = dimv[3] dimy = dimv[2] dimz = dimv[1] dimt = dimv[0] undimy = 'uny' unva = np.zeros((dimt,dimz,dimy-1,dimx), dtype=np.float) unva[...,0:dimy-1,:] = 0.5*(vunva[...,1:dimy,:]+vunva[...,0:dimy-1,:]) # Attributes of the variable Vvals = gen.variables_values('va') dictcompvars['unva'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} if UNwa_compute: print ' ' + main + ": un-staggering '" + vUNwa + "': as 0.5*(wa[0:dimz-1]+" + \ "wa[1:dimz])" vunwa = ncobj.variables[vUNwa][:] dimv = ncobj.variables[vUNwa].shape dimx = dimv[3] dimy = dimv[2] dimz = dimv[1] dimt = dimv[0] undimz = 'unz' unwa = np.zeros((dimt,dimz-1,dimy,dimx), dtype=np.float) unwa[...,0:dimz-1,:,:] = 0.5*(vunwa[...,1:dimz,:,:]+vunwa[...,0:dimz-1,:,:]) # Attributes of the variable Vvals = gen.variables_values('wa') dictcompvars['unwa'] = {'name': Vvals[0], 'standard_name': Vvals[1], \ 'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]} ### ## # # Going for the diagnostics ### ## # print ' ' + main + ' ...' varsadd = [] Varns = ncobj.variables.keys() Varns.sort() for idiag in range(Ndiags): print ' diagnostic:',diags[idiag] diagn = diags[idiag].split('|')[0] depvars = diags[idiag].split('|')[1].split('@') if not gen.searchInlist(NONcheckdepvars, diagn): if diags[idiag].split('|')[1].find('@') != -1: depvars = diags[idiag].split('|')[1].split('@') if depvars[0] == 'deaccum': diagn='deaccum' if depvars[0] == 'accum': diagn='accum' for depv in depvars: # Checking without extra arguments of a depending variable (':', separated) if depv.find(':') != -1: depv=depv.split(':')[0] if not ncobj.variables.has_key(depv) and not \ gen.searchInlist(NONcheckingvars, depv) and \ not gen.searchInlist(methods, depv) and not depvars[0] == 'deaccum'\ and not depvars[0] == 'accum' and not depv[0:2] == 'z=': print errormsg print ' ' + main + ": file '" + opts.ncfile + \ "' does not have variable '" + depv + "' !!" print ' available ones:', Varns quit(-1) else: depvars = diags[idiag].split('|')[1] if not ncobj.variables.has_key(depvars) and not \ gen.searchInlist(NONcheckingvars, depvars) and \ not gen.searchInlist(methods, depvars): print errormsg print ' ' + main + ": file '" + opts.ncfile + \ "' does not have variable '" + depvars + "' !!" print ' available ones:', Varns quit(-1) print "\n Computing '" + diagn + "' from: ", depvars, '...' # acraintot: accumulated total precipitation from WRF RAINC, RAINNC, RAINSH if diagn == 'ACRAINTOT': var0 = ncobj.variables[depvars[0]] var1 = ncobj.variables[depvars[1]] var2 = ncobj.variables[depvars[2]] diagout = var0[:] + var1[:] + var2[:] dnamesvar = var0.dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(dvnamesvar,nonvd): dvnamesvar.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'pracc', diagout, dnamesvar, dvnamesvar, newnc) # accum: acumulation of any variable as (Variable, time [as [tunits] # from/since ....], newvarname) elif diagn == 'accum': var0 = ncobj.variables[depvars[0]] var1 = ncobj.variables[depvars[1]] dnamesvar = var0.dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.compute_accum(var0,dnamesvar,dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) CFvarn = gen.variables_values(depvars[0])[0] # Removing the flux if depvars[1] == 'XTIME': dtimeunits = var1.getncattr('description') tunits = dtimeunits.split(' ')[0] else: dtimeunits = var1.getncattr('units') tunits = dtimeunits.split(' ')[0] dtime = (var1[1] - var1[0])*diag.timeunits_seconds(tunits) ncvar.insert_variable(ncobj, CFvarn + 'acc', diagout*dtime, diagoutd, diagoutvd, newnc) # cllmh with cldfra, pres elif diagn == 'cllmh': var0 = ncobj.variables[depvars[0]] if depvars[1] == 'WRFp': var1 = WRFp else: var01 = ncobj.variables[depvars[1]] if len(size(var1.shape)) < len(size(var0.shape)): var1 = np.brodcast_arrays(var01,var0)[0] else: var1 = var01 diagout, diagoutd, diagoutvd = diag.Forcompute_cllmh(var0,var1,dnames,dvnames) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'cll', diagout[0,:], diagoutd, diagoutvd, newnc) ncvar.insert_variable(ncobj, 'clm', diagout[1,:], diagoutd, diagoutvd, newnc) ncvar.insert_variable(ncobj, 'clh', diagout[2,:], diagoutd, diagoutvd, newnc) # clt with cldfra elif diagn == 'clt': var0 = ncobj.variables[depvars] diagout, diagoutd, diagoutvd = diag.Forcompute_clt(var0,dnames,dvnames) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'clt', diagout, diagoutd, diagoutvd, newnc) # convini (pr, time) elif diagn == 'convini': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] otime = ncobj.variables[depvars[1]] dnamesvar = ncobj.variables[depvars[0]].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.var_convini(var0, var1, dnames, dvnames) ncvar.insert_variable(ncobj, 'convini', diagout, diagoutd, diagoutvd, newnc, \ fill=gen.fillValueF) # Getting the right units ovar = newnc.variables['convini'] if gen.searchInlist(otime.ncattrs(), 'units'): tunits = otime.getncattr('units') ncvar.set_attribute(ovar, 'units', tunits) newnc.sync() # deaccum: deacumulation of any variable as (Variable, time [as [tunits] # from/since ....], newvarname) elif diagn == 'deaccum': var0 = ncobj.variables[depvars[0]] var1 = ncobj.variables[depvars[1]] dnamesvar = var0.dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.compute_deaccum(var0,dnamesvar,dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) # Transforming to a flux if depvars[1] == 'XTIME': dtimeunits = var1.getncattr('description') tunits = dtimeunits.split(' ')[0] else: dtimeunits = var1.getncattr('units') tunits = dtimeunits.split(' ')[0] dtime = (var1[1] - var1[0])*diag.timeunits_seconds(tunits) ncvar.insert_variable(ncobj, depvars[2], diagout/dtime, diagoutd, diagoutvd, \ newnc) # fog_K84: Computation of fog and visibility following Kunkel, (1984) as QCLOUD, QICE elif diagn == 'fog_K84': var0 = ncobj.variables[depvars[0]] var1 = ncobj.variables[depvars[1]] dnamesvar = list(var0.dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diag1, diag2, diagoutd, diagoutvd = diag.Forcompute_fog_K84(var0, var1, \ dnamesvar, dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'fog', diag1, diagoutd, diagoutvd, newnc) ncvar.insert_variable(ncobj, 'fogvisblty', diag2, diagoutd, diagoutvd, newnc) # fog_RUC: Computation of fog and visibility following Kunkel, (1984) as QVAPOR, # WRFt, WRFp or Q2, T2, PSFC elif diagn == 'fog_RUC': var0 = ncobj.variables[depvars[0]] print gen.infmsg if depvars[1] == 'WRFt': print ' ' + main + ": computing '" + diagn + "' using 3D variables !!" var1 = WRFt var2 = WRFp else: print ' ' + main + ": computing '" + diagn + "' using 2D variables !!" var1 = ncobj.variables[depvars[1]] var2 = ncobj.variables[depvars[2]] dnamesvar = list(var0.dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diag1, diag2, diagoutd, diagoutvd = diag.Forcompute_fog_RUC(var0, var1, var2,\ dnamesvar, dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'fog', diag1, diagoutd, diagoutvd, newnc) ncvar.insert_variable(ncobj, 'fogvisblty', diag2, diagoutd, diagoutvd, newnc) # fog_FRAML50: Computation of fog and visibility following Gultepe, I. and # J.A. Milbrandt, 2010 as QVAPOR, WRFt, WRFp or Q2, T2, PSFC elif diagn == 'fog_FRAML50': var0 = ncobj.variables[depvars[0]] print gen.infmsg if depvars[1] == 'WRFt': print ' ' + main + ": computing '" + diagn + "' using 3D variables !!" var1 = WRFt var2 = WRFp else: print ' ' + main + ": computing '" + diagn + "' using 2D variables !!" var1 = ncobj.variables[depvars[1]] var2 = ncobj.variables[depvars[2]] dnamesvar = list(var0.dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diag1, diag2, diagoutd, diagoutvd = diag.Forcompute_fog_FRAML50(var0, var1, \ var2, dnamesvar, dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'fog', diag1, diagoutd, diagoutvd, newnc) ncvar.insert_variable(ncobj, 'fogvisblty', diag2, diagoutd, diagoutvd, newnc) # front_R04: Computation of the presence of a front as defined by Rodrigues et al. # (2004) (tas, uas, vas, dxs, dys, 'params',[dtas],[dwss]) # Rodrigues et al. 2004: Climatologia de frentes frias no litoral de Santa Catarina, # Rev. Bras. Geof. v.22 n.2 Sao Paulo maio/ago. DOI: 10.1590/S0102-261X2004000200004 elif diagn == 'front_R04': var0 = ncobj.variables[depvars[0]] var1 = ncobj.variables[depvars[1]] var2 = ncobj.variables[depvars[2]] if depvars[3] == 'WRFdx' or depvars[3] == 'LONLATdx': var3 = WRFdx else: var3 = ncobj.variables[depvars[3]] if depvars[4] == 'WRFdy' or depvars[4] == 'LONLATdy': var4 = WRFdy else: var4 = ncobj.variables[depvars[4]] par1 = np.float(depvars[5].split(':')[1]) par2 = np.float(depvars[5].split(':')[2]) dnamesvar = list(var0.dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diag1, diag2, diag3, diag4, diagoutd, diagoutvd = diag.Forcompute_front_R04( \ var0, var1, var2, var3, var4, par1, par2, dnamesvar, dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'front', diag1, diagoutd, diagoutvd, newnc, \ fill=gen.fillValueF) ovar = newnc.variables['front'] ovar.setncattr('description', 'front defintion after a generalization of ' + \ 'Rodrigues et al. 2004, Rev. Bras. Geof. v.22 n.2, ' + \ 'doi: 10.1590/S0102-261X2004000200004') ovar.setncattr('details1', 'N-S wind gradient replaced by gradient of uas, '+\ 'vas') ovar.setncattr('dtas', par1) ovar.setncattr('dwss', par2) newnc.sync() ncvar.insert_variable(ncobj, 'dt1', diag2, diagoutd, diagoutvd, newnc, \ fill=gen.fillValueF) ovar = newnc.variables['dt1'] vu = var0.getncattr('units') ncvar.set_attribute(ovar, 'units', vu) attrv = ovar.getncattr('standard_name') ncvar.set_attribute(ovar, 'standard_name', attrv + depvars[0]) attrv = ovar.getncattr('long_name') ncvar.set_attribute(ovar, 'long_name', attrv + ' of ' + depvars[0]) ncobj.sync() ncvar.insert_variable(ncobj, 'gradh', diag3, diagoutd, diagoutvd, newnc, \ fill=gen.fillValueF) ovar = newnc.variables['gradh'] vu = var1.getncattr('units') ncvar.set_attribute(ovar, 'units', vu+'m-1') attrv = ovar.getncattr('standard_name') ncvar.set_attribute(ovar, 'standard_name', attrv + 'wss') attrv = ovar.getncattr('long_name') ncvar.set_attribute(ovar, 'long_name', attrv + ' of wss') ncobj.sync() ncvar.insert_variable(ncobj, 'dt2', diag4, diagoutd, diagoutvd, newnc, \ fill=gen.fillValueF) ovar = newnc.variables['dt2'] vu = var0.getncattr('units') ncvar.set_attribute(ovar, 'units', vu) attrv = ovar.getncattr('standard_name') ncvar.set_attribute(ovar, 'standard_name', attrv + depvars[0]) attrv = ovar.getncattr('long_name') ncvar.set_attribute(ovar, 'long_name', attrv + ' of ' + depvars[0]) ncobj.sync() # frontogenesis (theta, ua, va, wa, press, dxs, dys, dzs, time) elif diagn == 'frontogenesis': if depvars[0] == 'WRFt': var0 = WRFt else: var0 = ncobj.variables[depvars[0]][:] dx = var0.shape[3] dy = var0.shape[2] dz = var0.shape[1] dt = var0.shape[0] if depvars[1][0:4] == 'UNua': var1 = unua else: var1 = ncobj.variables[depvars[1]][:] if depvars[2][0:4] == 'UNva': var2 = unva else: var2 = ncobj.variables[depvars[2]][:] if depvars[3][0:4] == 'UNwa': var3 = unwa else: var3 = ncobj.variables[depvars[3]][:] if depvars[4] == 'WRFp': var4 = WRFp else: var4 = ncobj.variables[depvars[4]][:] if depvars[5] == 'WRFdx': var5 = WRFdx else: var5 = ncobj.variables[depvars[5]][:] if depvars[6] == 'WRFdy': var6 = WRFdy else: var6 = ncobj.variables[depvars[6]][:] if depvars[7] == 'WRFdz': var7 = WRFdz[0,0:dz,0,0] else: var7 = ncobj.variables[depvars[7]][:] if depvars[8] == 'WRFtime': var08 = WRFtime else: var08 = ncobj.variables[depvars[8]][:] # Assuming monotonic time-axis... var8 = var08[1] - var08[0] if len(var4.shape) == 1: print ' ' + diagn + ': rank 1 press !!' print ' spreading over 4D !!' var40 = var4 + 0. var4 = np.zeros((dt,dz,dy,dx), dtype=np.float) for it in range(dt): for iy in range(dy): for ix in range(dx): var4[it,:,iy,ix] = var40 # Unifying... #var5 = var5/var5 #var6 = var6/var6 #var7 = var7/var7 diag1, diag2, diag3, diag4, diag5, diag6, diag7, diag8, diag9, diag10, \ diagoutd, diagoutvd = diag.Forcompute_frontogenesis(var0, var1, var2, var3,\ var4, var5, var6, var7, var8, dnames, dvnames) # Removing inf, Nan, .... diag1 = ma.masked_invalid(diag1) diag2 = ma.masked_invalid(diag2) diag3 = ma.masked_invalid(diag3) diag4 = ma.masked_invalid(diag4) diag5 = ma.masked_invalid(diag5) diag6 = ma.masked_invalid(diag6) diag7 = ma.masked_invalid(diag7) diag8 = ma.masked_invalid(diag8) diag9 = ma.masked_invalid(diag9) diag10 = ma.masked_invalid(diag10) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'diabh', diag1, diagoutd, diagoutvd, newnc, \ gen.fillValueF) newnc.renameVariable('diabh', 'xdiabh') newnc.sync() ovar = newnc.variables['xdiabh'] stdn = ovar.getncattr('standard_name') ovar.setncattr('standard_name', 'x'+stdn) lngn = ovar.getncattr('long_name') ovar.setncattr('long_name', 'x component of theta '+stdn) uts = ovar.getncattr('units') ovar.setncattr('units', 'Km-1s-1') newnc.sync() ncvar.insert_variable(ncobj, 'diabh', diag2, diagoutd, diagoutvd, newnc, \ gen.fillValueF) newnc.renameVariable('diabh', 'ydiabh') newnc.sync() ovar = newnc.variables['ydiabh'] stdn = ovar.getncattr('standard_name') ovar.setncattr('standard_name', 'y'+stdn) lngn = ovar.getncattr('long_name') ovar.setncattr('long_name', 'y component of theta '+stdn) uts = ovar.getncattr('units') ovar.setncattr('units', 'Km-1s-1') newnc.sync() ncvar.insert_variable(ncobj, 'diabh', diag3, diagoutd, diagoutvd, newnc, \ gen.fillValueF) newnc.renameVariable('diabh', 'zdiabh') newnc.sync() ovar = newnc.variables['zdiabh'] stdn = ovar.getncattr('standard_name') ovar.setncattr('standard_name', 'z'+stdn) lngn = ovar.getncattr('long_name') ovar.setncattr('long_name', 'z component of theta '+stdn) uts = ovar.getncattr('units') ovar.setncattr('units', 'Km-1s-1') newnc.sync() ncvar.insert_variable(ncobj, 'def', diag4, diagoutd, diagoutvd, newnc, \ gen.fillValueF) newnc.renameVariable('def', 'xdef') newnc.sync() ovar = newnc.variables['xdef'] stdn = ovar.getncattr('standard_name') ovar.setncattr('standard_name', 'thetax'+stdn) lngn = ovar.getncattr('long_name') ovar.setncattr('long_name', 'x component of theta '+stdn) uts = ovar.getncattr('units') ovar.setncattr('units', 'Km-1s-1') newnc.sync() ncvar.insert_variable(ncobj, 'def', diag5, diagoutd, diagoutvd, newnc, \ gen.fillValueF) newnc.renameVariable('def', 'ydef') newnc.sync() ovar = newnc.variables['ydef'] stdn = ovar.getncattr('standard_name') ovar.setncattr('standard_name', 'thetay'+stdn) lngn = ovar.getncattr('long_name') ovar.setncattr('long_name', 'y component of theta '+stdn) uts = ovar.getncattr('units') ovar.setncattr('units', 'Km-1s-1') newnc.sync() ncvar.insert_variable(ncobj, 'def', diag6, diagoutd, diagoutvd, newnc, \ gen.fillValueF) newnc.renameVariable('def', 'zdef') newnc.sync() ovar = newnc.variables['zdef'] stdn = ovar.getncattr('standard_name') ovar.setncattr('standard_name', 'thetaz'+stdn) lngn = ovar.getncattr('long_name') ovar.setncattr('long_name', 'z component of theta '+stdn) uts = ovar.getncattr('units') ovar.setncattr('units', 'Km-1s-1') newnc.sync() ncvar.insert_variable(ncobj, 'tilt', diag7, diagoutd, diagoutvd, newnc, \ gen.fillValueF) newnc.renameVariable('tilt', 'xtilt') newnc.sync() ovar = newnc.variables['xtilt'] stdn = ovar.getncattr('standard_name') ovar.setncattr('standard_name', 'thetax'+stdn) lngn = ovar.getncattr('long_name') ovar.setncattr('long_name', 'x component of theta '+stdn) uts = ovar.getncattr('units') ovar.setncattr('units', 'Km-1s-1') newnc.sync() ncvar.insert_variable(ncobj, 'tilt', diag8, diagoutd, diagoutvd, newnc, \ gen.fillValueF) newnc.renameVariable('tilt', 'ytilt') newnc.sync() ovar = newnc.variables['ytilt'] stdn = ovar.getncattr('standard_name') ovar.setncattr('standard_name', 'thetay'+stdn) lngn = ovar.getncattr('long_name') ovar.setncattr('long_name', 'y component of theta '+stdn) uts = ovar.getncattr('units') ovar.setncattr('units', 'Km-1s-1') newnc.sync() ncvar.insert_variable(ncobj, 'div', diag9, diagoutd, diagoutvd, newnc, \ gen.fillValueF) newnc.renameVariable('div', 'zdiv') newnc.sync() ovar = newnc.variables['zdiv'] stdn = ovar.getncattr('standard_name') ovar.setncattr('standard_name', 'thetaz'+stdn) lngn = ovar.getncattr('long_name') ovar.setncattr('long_name', 'x component of theta '+stdn) uts = ovar.getncattr('units') ovar.setncattr('units', 'Km-1s-1') newnc.sync() newdim = newnc.createDimension('e', 3) newvar = newnc.createVariable('e', 'i', ('e')) newvar[:] = range(3) ncvar.basicvardef(newvar, 'component', 'axis component: x, y, z', '-') newnc.sync() diagoutd = ['e'] + diagoutd diagoutvd = ['e'] + diagoutvd ncvar.insert_variable(ncobj, 'frontogenesis', diag10, diagoutd, diagoutvd, \ newnc, gen.fillValueF) newnc.sync() # gradient2Dh: 1st order Horizontal 2D-gradient of any variable [variable], [distx], # [disty] elif diagn == 'gradient2Dh': var0 = ncobj.variables[depvars[0]] if depvars[1] == 'WRFdx': var1 = WRFdx else: var1 = ncobj.variables[depvars[1]][0,:,:] if depvars[2] == 'WRFdy': var2 = WRFdy else: var2 = ncobj.variables[depvars[2]][0,:,:] dnamesvar = var0.dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.Forcompute_gradient2Dh(var0, var1, var2, \ dnamesvar, dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) CFvarn = gen.variables_values(depvars[0])[0] ncvar.insert_variable(ncobj, depvars[0], diagout, diagoutd, diagoutvd, \ newnc, gen.fillValueF) newnc.sync() CFvar = gen.variables_values(depvars[0]) newnc.renameVariable(CFvar[0], CFvar[0]+'grad2dh') newnc.sync() ovar = newnc.variables[CFvar[0]+'grad2dh' varu = ovar.units ovar.setncattr('units', varu+'ds-1') newnc.sync() # LMDZrh (pres, t, r) elif diagn == 'LMDZrh': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] var2 = ncobj.variables[depvars[2]][:] diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnames,dvnames) ncvar.insert_variable(ncobj, 'hur', diagout, diagoutd, diagoutvd, newnc) # LMDZrhs (psol, t2m, q2m) elif diagn == 'LMDZrhs': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] var2 = ncobj.variables[depvars[2]][:] dnamesvar = ncobj.variables[depvars[0]].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar) ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc) # range_faces: LON, LAT, HGT, WRFdxdy, 'face:['WE'/'SN']:[dsfilt]:[dsnewrange]:[hvalleyrange]' elif diagn == 'range_faces': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] var2 = ncobj.variables[depvars[2]][:] face = depvars[4].split(':')[1] dsfilt = np.float(depvars[4].split(':')[2]) dsnewrange = np.float(depvars[4].split(':')[3]) hvalleyrange = np.float(depvars[4].split(':')[4]) dnamesvar = list(ncobj.variables[depvars[2]].dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) lon, lat = gen.lonlat2D(var0, var1) if len(var2.shape) == 3: print warnmsg print ' ' + diagn + ": shapping to 2D variable '" + depvars[2] + "' !!" hgt = var2[0,:,:] dnamesvar.pop(0) dvnamesvar.pop(0) else: hgt = var2[:] orogmax, ptorogmax, dhgt, peaks, valleys, origfaces, diagout, diagoutd, \ diagoutvd, rng, rngorogmax, ptrngorogmax= diag.Forcompute_range_faces(lon, \ lat, hgt, WRFdx, WRFdy, WRFds, face, dsfilt, dsnewrange, hvalleyrange, \ dnamesvar, dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'dx', WRFdx, diagoutd, diagoutvd, newnc) ncvar.insert_variable(ncobj, 'dy', WRFdy, diagoutd, diagoutvd, newnc) ncvar.insert_variable(ncobj, 'ds', WRFds, diagoutd, diagoutvd, newnc) # adding variables to output file if face == 'WE': axis = 'lon' elif face == 'SN': axis = 'lat' ncvar.insert_variable(ncobj, 'range', rng, diagoutd, diagoutvd, newnc, \ fill=gen.fillValueI) ovar = newnc.variables['range'] ncvar.set_attribute(ovar, 'deriv', axis) ncvar.insert_variable(ncobj, 'orogmax', rngorogmax, diagoutd, diagoutvd, \ newnc, fill=gen.fillValueF) newnc.renameVariable('orogmax', 'rangeorogmax') ovar = newnc.variables['rangeorogmax'] ncvar.set_attribute(ovar, 'deriv', axis) stdn = ovar.standard_name ncvar.set_attribute(ovar, 'standard_name', 'range_' + stdn) Ln = ovar.long_name ncvar.set_attribute(ovar, 'long_name', 'range ' + stdn) ncvar.insert_variable(ncobj, 'ptorogmax', ptrngorogmax, diagoutd, diagoutvd, \ newnc) newnc.renameVariable('ptorogmax', 'rangeptorogmax') ovar = newnc.variables['rangeptorogmax'] ncvar.set_attribute(ovar, 'deriv', axis) stdn = ovar.standard_name ncvar.set_attribute(ovar, 'standard_name', 'range_' + stdn) Ln = ovar.long_name ncvar.set_attribute(ovar, 'long_name', 'range ' + stdn) ncvar.insert_variable(ncobj, 'orogmax', orogmax, diagoutd, diagoutvd, \ newnc) ovar = newnc.variables['orogmax'] ncvar.set_attribute(ovar, 'deriv', axis) ncvar.insert_variable(ncobj, 'ptorogmax', ptorogmax, diagoutd, diagoutvd, \ newnc) ovar = newnc.variables['ptorogmax'] ncvar.set_attribute(ovar, 'deriv', axis) ncvar.insert_variable(ncobj, 'orogderiv', dhgt, diagoutd, diagoutvd, newnc) ovar = newnc.variables['orogderiv'] ncvar.set_attribute(ovar, 'deriv', axis) ncvar.insert_variable(ncobj, 'peak', peaks, diagoutd, diagoutvd, newnc) ncvar.insert_variable(ncobj, 'valley', valleys, diagoutd, diagoutvd, newnc) ncvar.insert_variable(ncobj, 'rangefaces', diagout, diagoutd, diagoutvd, \ newnc) newnc.renameVariable('rangefaces', 'rangefacesfilt') ovar = newnc.variables['rangefacesfilt'] ncvar.set_attribute(ovar, 'face', face) ncvar.set_attributek(ovar, 'dist_filter', dsfilt, 'F') ncvar.insert_variable(ncobj, 'rangefaces', origfaces, diagoutd, diagoutvd, \ newnc, fill=gen.fillValueI) ovar = newnc.variables['rangefaces'] ncvar.set_attribute(ovar, 'face', face) ncvar.set_attributek(ovar, 'dist_newrange', dsnewrange, 'F') ncvar.set_attributek(ovar, 'h_valley_newrange', hvalleyrange, 'F') # cell_bnds: grid cell bounds from lon, lat from a reglar lon/lat projection as # intersection of their related parallels and meridians elif diagn == 'reglonlatbnds': var00 = ncobj.variables[depvars[0]][:] var01 = ncobj.variables[depvars[1]][:] var0, var1 = gen.lonlat2D(var00,var01) dnamesvar = [] dnamesvar.append('bnds') if (len(var00.shape) == 3): dnamesvar.append(ncobj.variables[depvars[0]].dimensions[1]) dnamesvar.append(ncobj.variables[depvars[0]].dimensions[2]) elif (len(var00.shape) == 2): dnamesvar.append(ncobj.variables[depvars[0]].dimensions[0]) dnamesvar.append(ncobj.variables[depvars[0]].dimensions[1]) elif (len(var00.shape) == 1): dnamesvar.append(ncobj.variables[depvars[0]].dimensions[0]) dnamesvar.append(ncobj.variables[depvars[1]].dimensions[0]) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) cellbndsx, cellbndsy, diagoutd, diagoutvd = diag.Forcompute_cellbndsreg(var0,\ var1, dnamesvar, dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) # creation of bounds dimension newdim = newnc.createDimension('bnds', 4) ncvar.insert_variable(ncobj, 'lon_bnds', cellbndsx, diagoutd, diagoutvd, newnc) ncvar.insert_variable(ncobj, 'lat_bnds', cellbndsy, diagoutd, diagoutvd, newnc) # tws: Bet Wulb temperature following Stull 2011 (tas, hurs) elif diagn == 'tws': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] dnamesvar = list(ncobj.variables[depvars[0]].dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagoutd = dnames diagoutvd = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout = diag.var_tws_S11(var0,var1) ncvar.insert_variable(ncobj, 'tws', diagout, diagoutd, diagoutvd, newnc) # cell_bnds: grid cell bounds from XLONG_U, XLAT_U, XLONG_V, XLAT_V as intersection # of their related parallels and meridians elif diagn == 'WRFbnds': var0 = ncobj.variables[depvars[0]][0,:,:] var1 = ncobj.variables[depvars[1]][0,:,:] var2 = ncobj.variables[depvars[2]][0,:,:] var3 = ncobj.variables[depvars[3]][0,:,:] dnamesvar = [] dnamesvar.append('bnds') dnamesvar.append(ncobj.variables[depvars[0]].dimensions[1]) dnamesvar.append(ncobj.variables[depvars[2]].dimensions[2]) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) cellbndsx, cellbndsy, diagoutd, diagoutvd = diag.Forcompute_cellbnds(var0, \ var1, var2, var3, dnamesvar, dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) # creation of bounds dimension newdim = newnc.createDimension('bnds', 4) ncvar.insert_variable(ncobj, 'lon_bnds', cellbndsx, diagoutd, diagoutvd, newnc) newnc.sync() ncvar.insert_variable(ncobj, 'lat_bnds', cellbndsy, diagoutd, diagoutvd, newnc) newnc.sync() if newnc.variables.has_key('XLONG'): ovar = newnc.variables['XLONG'] ovar.setncattr('bounds', 'lon_bnds lat_bnds') ovar = newnc.variables['XLAT'] ovar.setncattr('bounds', 'lon_bnds lat_bnds') elif newnc.variables.has_key('XLONG_M'): ovar = newnc.variables['XLONG_M'] ovar.setncattr('bounds', 'lon_bnds lat_bnds') ovar = newnc.variables['XLAT_M'] ovar.setncattr('bounds', 'lon_bnds lat_bnds') else: print errormsg print ' ' + fname + ": error computing diagnostic '" + diagn + "' !!" print " neither: 'XLONG', 'XLONG_M' have been found" quit(-1) # mrso: total soil moisture SMOIS, DZS elif diagn == 'WRFmrso': var0 = ncobj.variables[depvars[0]][:] var10 = ncobj.variables[depvars[1]][:] dnamesvar = list(ncobj.variables[depvars[0]].dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) var1 = var0.copy()*0. var2 = var0.copy()*0.+1. # Must be a better way.... for j in range(var0.shape[2]): for i in range(var0.shape[3]): var1[:,:,j,i] = var10 diagout, diagoutd, diagoutvd = diag.Forcompute_zint(var0, var1, var2, \ dnamesvar, dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'mrso', diagout, diagoutd, diagoutvd, newnc) # mrsos: First layer soil moisture SMOIS, DZS elif diagn == 'WRFmrsos': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] diagoutd = list(ncobj.variables[depvars[0]].dimensions) diagoutvd = ncvar.var_dim_dimv(diagoutd,dnames,dvnames) diagoutd.pop(1) diagoutvd.pop(1) diagout= np.zeros((var0.shape[0],var0.shape[2],var0.shape[3]), dtype=np.float) # Must be a better way.... for j in range(var0.shape[2]): for i in range(var0.shape[3]): diagout[:,j,i] = var0[:,0,j,i]*var1[:,0] # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'mrsos', diagout, diagoutd, diagoutvd, newnc) # mslp: mean sea level pressure (pres, psfc, terrain, temp, qv) elif diagn == 'mslp' or diagn == 'WRFmslp': var1 = ncobj.variables[depvars[1]][:] var2 = ncobj.variables[depvars[2]][:] var4 = ncobj.variables[depvars[4]][:] if diagn == 'WRFmslp': var0 = WRFp var3 = WRFt dnamesvar = ncobj.variables['P'].dimensions else: var0 = ncobj.variables[depvars[0]][:] var3 = ncobj.variables[depvars[3]][:] dnamesvar = ncobj.variables[depvars[0]].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.compute_mslp(var0, var1, var2, var3, var4, \ dnamesvar, dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc) # WRFtws: Bet Wulb temperature following Stull 2011 (PSFC, T2, Q2) elif diagn == 'WRFtws': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] var2 = ncobj.variables[depvars[2]][:] dnamesvar = list(ncobj.variables[depvars[2]].dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) hurs, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar) diagout = diag.var_tws_S11(var1, hurs) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'tws', diagout, diagoutd, diagoutvd, newnc) # OMEGAw (omega, p, t) from NCL formulation (https://www.ncl.ucar.edu/Document/Functions/Contributed/omega_to_w.shtml) elif diagn == 'OMEGAw': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] var2 = ncobj.variables[depvars[2]][:] dnamesvar = ncobj.variables[depvars[0]].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.compute_OMEGAw(var0,var1,var2,dnamesvar,dvnamesvar) ncvar.insert_variable(ncobj, 'wa', diagout, diagoutd, diagoutvd, newnc) # raintot: instantaneous total precipitation from WRF as (RAINC + RAINC + RAINSH) / dTime elif diagn == 'RAINTOT': var0 = ncobj.variables[depvars[0]] var1 = ncobj.variables[depvars[1]] var2 = ncobj.variables[depvars[2]] if depvars[3] != 'WRFtime': var3 = ncobj.variables[depvars[3]] else: var3 = np.arange(var0.shape[0], dtype=int) var = var0[:] + var1[:] + var2[:] dnamesvar = var0.dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.compute_deaccum(var,dnamesvar,dvnamesvar) # Transforming to a flux if var3.shape[0] > 1: if depvars[3] != 'WRFtime': dtimeunits = var3.getncattr('units') tunits = dtimeunits.split(' ')[0] dtime = (var3[1] - var3[0])*diag.timeunits_seconds(tunits) else: var3 = ncobj.variables['Times'] time1 = var3[0,:] time2 = var3[1,:] tmf1 = '' tmf2 = '' for ic in range(len(time1)): tmf1 = tmf1 + time1[ic] tmf2 = tmf2 + time2[ic] dtdate1 = dtime.datetime.strptime(tmf1,"%Y-%m-%d_%H:%M:%S") dtdate2 = dtime.datetime.strptime(tmf2,"%Y-%m-%d_%H:%M:%S") diffdate12 = dtdate2 - dtdate1 dtime = diffdate12.total_seconds() print 'dtime:',dtime else: print warnmsg print ' ' + main + ": only 1 time-step for '" + diag + "' !!" print ' leaving a zero value!' diagout = var0[:]*0. dtime=1. # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'pr', diagout/dtime, diagoutd, diagoutvd, newnc) # timemax ([varname], time). When a given variable [varname] got its maximum elif diagn == 'timemax': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] otime = ncobj.variables[depvars[1]] dnamesvar = ncobj.variables[depvars[0]].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.var_timemax(var0, var1, dnames, \ dvnames) ncvar.insert_variable(ncobj, 'timemax', diagout, diagoutd, diagoutvd, newnc, \ fill=gen.fillValueF) # Getting the right units ovar = newnc.variables['timemax'] if gen.searchInlist(otime.ncattrs(), 'units'): tunits = otime.getncattr('units') ncvar.set_attribute(ovar, 'units', tunits) newnc.sync() ncvar.set_attribute(ovar, 'variable', depvars[0]) # timeoverthres ([varname], time, [value], [CFvarn]). When a given variable [varname] # overpass a given [value]. Being [CFvarn] the name of the diagnostics in # variables_values.dat elif diagn == 'timeoverthres': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] var2 = np.float(depvars[2]) var3 = depvars[3] otime = ncobj.variables[depvars[1]] dnamesvar = ncobj.variables[depvars[0]].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.var_timeoverthres(var0, var1, var2, \ dnames, dvnames) ncvar.insert_variable(ncobj, var3, diagout, diagoutd, diagoutvd, newnc, \ fill=gen.fillValueF) # Getting the right units ovar = newnc.variables[var3] if gen.searchInlist(otime.ncattrs(), 'units'): tunits = otime.getncattr('units') ncvar.set_attribute(ovar, 'units', tunits) newnc.sync() ncvar.set_attribute(ovar, 'overpassed_threshold', var2) # rhs (psfc, t, q) from TimeSeries files elif diagn == 'TSrhs': p0=100000. var0 = ncobj.variables[depvars[0]][:] var1 = (ncobj.variables[depvars[1]][:])*(var0/p0)**(2./7.) var2 = ncobj.variables[depvars[2]][:] dnamesvar = ncobj.variables[depvars[0]].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar) ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc) # rhs (psfc, t, q) from tas, tds elif diagn == 'rhs_tas_tds': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] dnamesvar = ncobj.variables[depvars[0]].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.var_hur_tas_tds(var0,var1,dnamesvar, \ dvnamesvar) ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc) # slw: total soil liquid water SH2O, DZS elif diagn == 'WRFslw': var0 = ncobj.variables[depvars[0]][:] var10 = ncobj.variables[depvars[1]][:] dnamesvar = list(ncobj.variables[depvars[0]].dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) var1 = var0.copy()*0. var2 = var0.copy()*0.+1. # Must be a better way.... for j in range(var0.shape[2]): for i in range(var0.shape[3]): var1[:,:,j,i] = var10 diagout, diagoutd, diagoutvd = diag.Forcompute_zint(var0, var1, var2, \ dnamesvar, dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'slw', diagout, diagoutd, diagoutvd, newnc) # td (psfc, t, q) from TimeSeries files elif diagn == 'TStd' or diagn == 'td': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] - 273.15 var2 = ncobj.variables[depvars[2]][:] dnamesvar = ncobj.variables[depvars[0]].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar) ncvar.insert_variable(ncobj, 'tdas', diagout, diagoutd, diagoutvd, newnc) # td (psfc, t, q) from TimeSeries files elif diagn == 'TStdC' or diagn == 'tdC': var0 = ncobj.variables[depvars[0]][:] # Temperature is already in degrees Celsius var1 = ncobj.variables[depvars[1]][:] var2 = ncobj.variables[depvars[2]][:] dnamesvar = ncobj.variables[depvars[0]].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'tdas', diagout, diagoutd, diagoutvd, newnc) # wds (u, v) elif diagn == 'TSwds' or diagn == 'wds' : var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] dnamesvar = ncobj.variables[depvars[0]].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.compute_wds(var0,var1,dnamesvar,dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'wds', diagout, diagoutd, diagoutvd, newnc) # wss (u, v) elif diagn == 'TSwss' or diagn == 'wss': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] dnamesvar = ncobj.variables[depvars[0]].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.compute_wss(var0,var1,dnamesvar,dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'wss', diagout, diagoutd, diagoutvd, newnc) # turbulence (var) elif diagn == 'turbulence': var0 = ncobj.variables[depvars][:] dnamesvar = list(ncobj.variables[depvars].dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.compute_turbulence(var0,dnamesvar,dvnamesvar) valsvar = gen.variables_values(depvars) newvarn = depvars + 'turb' ncvar.insert_variable(ncobj, newvarn, diagout, diagoutd, diagoutvd, newnc) varobj = newnc.variables[newvarn] attrv = varobj.long_name attr = varobj.delncattr('long_name') newattr = ncvar.set_attribute(varobj, 'long_name', attrv + \ " Taylor decomposition turbulence term") # ua va from ws wd (deg) elif diagn == 'uavaFROMwswd': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] ua = var0*np.cos(var1*np.pi/180.) va = var0*np.sin(var1*np.pi/180.) dnamesvar = ncobj.variables[depvars[0]].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) ncvar.insert_variable(ncobj, 'ua', ua, dnamesvar, dvnamesvar, newnc) ncvar.insert_variable(ncobj, 'va', va, dnamesvar, dvnamesvar, newnc) # ua va from obs ws wd (deg) elif diagn == 'uavaFROMobswswd': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] ua = var0*np.cos((var1+180.)*np.pi/180.) va = var0*np.sin((var1+180.)*np.pi/180.) dnamesvar = ncobj.variables[depvars[0]].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) ncvar.insert_variable(ncobj, 'ua', ua, dnamesvar, dvnamesvar, newnc) ncvar.insert_variable(ncobj, 'va', va, dnamesvar, dvnamesvar, newnc) # WRFbils fom WRF as HFX + LH elif diagn == 'WRFbils': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] diagout = var0 + var1 dnamesvar = list(ncobj.variables[depvars[0]].dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) ncvar.insert_variable(ncobj, 'bils', diagout, dnamesvar, dvnamesvar, newnc) # WRFcape_afwa CAPE, CIN, ZLFC, PLFC, LI following WRF 'phys/module_diaf_afwa.F' # methodology as WRFt, WRFrh, WRFp, WRFgeop, HGT elif diagn == 'WRFcape_afwa': var0 = WRFt var1 = WRFrh var2 = WRFp dz = WRFgeop.shape[1] # de-staggering var3 = 0.5*(WRFgeop[:,0:dz-1,:,:]+WRFgeop[:,1:dz,:,:])/9.8 var4 = ncobj.variables[depvars[4]][0,:,:] dnamesvar = list(ncobj.variables['T'].dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout = np.zeros(var0.shape, dtype=np.float) diagout1, diagout2, diagout3, diagout4, diagout5, diagoutd, diagoutvd = \ diag.Forcompute_cape_afwa(var0, var1, var2, var3, var4, 3, dnamesvar, \ dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'cape', diagout1, diagoutd, diagoutvd, newnc) ncvar.insert_variable(ncobj, 'cin', diagout2, diagoutd, diagoutvd, newnc) ncvar.insert_variable(ncobj, 'zlfc', diagout3, diagoutd, diagoutvd, newnc) ncvar.insert_variable(ncobj, 'plfc', diagout4, diagoutd, diagoutvd, newnc) ncvar.insert_variable(ncobj, 'li', diagout5, diagoutd, diagoutvd, newnc) # WRFclivi WRF water vapour path WRFdens, QICE, QGRAUPEL, QHAIL elif diagn == 'WRFclivi': var0 = WRFdens qtot = ncobj.variables[depvars[1]] qtotv = qtot[:] Nspecies = len(depvars) - 2 for iv in range(Nspecies): if ncobj.variables.has_key(depvars[iv+2]): var1 = ncobj.variables[depvars[iv+2]][:] qtotv = qtotv + var1 dnamesvar = list(qtot.dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.compute_clivi(var0, qtotv, dnamesvar,dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'clivi', diagout, diagoutd, diagoutvd, newnc) # WRFclwvi WRF water cloud-condensed path WRFdens, QCLOUD, QICE, QGRAUPEL, QHAIL elif diagn == 'WRFclwvi': var0 = WRFdens qtot = ncobj.variables[depvars[1]] qtotv = ncobj.variables[depvars[1]] Nspecies = len(depvars) - 2 for iv in range(Nspecies): if ncobj.variables.has_key(depvars[iv+2]): var1 = ncobj.variables[depvars[iv+2]] qtotv = qtotv + var1[:] dnamesvar = list(qtot.dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.compute_clwvl(var0, qtotv, dnamesvar,dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'clwvi', diagout, diagoutd, diagoutvd, newnc) # WRF_denszint WRF vertical integration as WRFdens, sum(Q[water species1], ..., Q[water speciesN]), varn=[varN] elif diagn == 'WRF_denszint': var0 = WRFdens varn = depvars[1].split('=')[1] qtot = ncobj.variables[depvars[2]] qtotv = ncobj.variables[depvars[2]] Nspecies = len(depvars) - 2 for iv in range(Nspecies): if ncobj.variables.has_key(depvars[iv+2]): var1 = ncobj.variables[depvars[iv+2]] qtotv = qtotv + var1[:] dnamesvar = list(qtot.dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.compute_clwvl(var0, qtotv, dnamesvar,dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, varn, diagout, diagoutd, diagoutvd, newnc) # WRFgeop geopotential from WRF as PH + PHB elif diagn == 'WRFgeop': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] # de-staggering geopotential diagout0 = var0 + var1 dt = diagout0.shape[0] dz = diagout0.shape[1] dy = diagout0.shape[2] dx = diagout0.shape[3] diagout = np.zeros((dt,dz-1,dy,dx), dtype=np.float) diagout = 0.5*(diagout0[:,1:dz,:,:]+diagout0[:,0:dz-1,:,:]) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'zg', diagout, dnames, diagoutvd, newnc) # WRFpotevap_orPM potential evapotranspiration following Penman-Monteith formulation # implemented in ORCHIDEE (in src_sechiba/enerbil.f90) as: WRFdens, UST, U10, V10, T2, PSFC, QVAPOR elif diagn == 'WRFpotevap_orPM': var0 = WRFdens[:,0,:,:] var1 = ncobj.variables[depvars[1]][:] var2 = ncobj.variables[depvars[2]][:] var3 = ncobj.variables[depvars[3]][:] var4 = ncobj.variables[depvars[4]][:] var5 = ncobj.variables[depvars[5]][:] var6 = ncobj.variables[depvars[6]][:,0,:,:] dnamesvar = list(ncobj.variables[depvars[1]].dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout = np.zeros(var1.shape, dtype=np.float) diagout, diagoutd, diagoutvd = diag.Forcompute_potevap_orPM(var0, var1, var2,\ var3, var4, var5, var6, dnamesvar, dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'evspsblpot', diagout, diagoutd, diagoutvd, newnc) # WRFmslp_ptarget sea-level pressure following ECMWF method as PSFC, HGT, WRFt, WRFp, ZNU, ZNW elif diagn == 'WRFpsl_ecmwf': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][0,:,:] var2 = WRFt[:,0,:,:] var4 = WRFp[:,0,:,:] var5 = ncobj.variables[depvars[4]][0,:] var6 = ncobj.variables[depvars[5]][0,:] # This is quite too appriximate!! passing pressure at half-levels to 2nd full # level, using eta values at full (ZNW) and half (ZNU) mass levels var3 = WRFp[:,0,:,:] + (var6[1] - var5[0])*(WRFp[:,1,:,:] - WRFp[:,0,:,:])/ \ (var5[1]-var5[0]) dnamesvar = list(ncobj.variables[depvars[0]].dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout = np.zeros(var0.shape, dtype=np.float) diagout, diagoutd, diagoutvd = diag.Forcompute_psl_ecmwf(var0, var1, var2, \ var3, var4, dnamesvar, dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc) # WRFmslp_ptarget sea-level pressure following ptarget method as WRFp, PSFC, WRFt, HGT, QVAPOR elif diagn == 'WRFpsl_ptarget': var0 = WRFp var1 = ncobj.variables[depvars[1]][:] var2 = WRFt var3 = ncobj.variables[depvars[3]][0,:,:] var4 = ncobj.variables[depvars[4]][:] dnamesvar = list(ncobj.variables[depvars[4]].dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout = np.zeros(var0.shape, dtype=np.float) diagout, diagoutd, diagoutvd = diag.Forcompute_psl_ptarget(var0, var1, var2, \ var3, var4, 700000., dnamesvar, dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc) # WRFp pressure from WRF as P + PB elif diagn == 'WRFp': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] diagout = var0 + var1 diagoutd = list(ncobj.variables[depvars[0]].dimensions) diagoutvd = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'pres', diagout, diagoutd, diagoutvd, newnc) # WRFpos elif diagn == 'WRFpos': dnamesvar = ncobj.variables['MAPFAC_M'].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) ncvar.insert_variable(ncobj, 'WRFpos', WRFpos, dnamesvar, dvnamesvar, newnc) # WRFprw WRF water vapour path WRFdens, QVAPOR elif diagn == 'WRFprw': var0 = WRFdens var1 = ncobj.variables[depvars[1]] dnamesvar = list(var1.dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.compute_prw(var0, var1, dnamesvar,dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'prw', diagout, diagoutd, diagoutvd, newnc) # WRFrh (P, T, QVAPOR) elif diagn == 'WRFrh': dnamesvar = list(ncobj.variables[depvars[2]].dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) ncvar.insert_variable(ncobj, 'hur', WRFrh, dnames, dvnames, newnc) # WRFrhs (PSFC, T2, Q2) elif diagn == 'WRFrhs': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] var2 = ncobj.variables[depvars[2]][:] dnamesvar = list(ncobj.variables[depvars[2]].dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar) ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc) # rvors (u10, v10, WRFpos) elif diagn == 'WRFrvors': var0 = ncobj.variables[depvars[0]] var1 = ncobj.variables[depvars[1]] diagout = rotational_z(var0, var1, distx) dnamesvar = ncobj.variables[depvars[0]].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) ncvar.insert_variable(ncobj, 'rvors', diagout, dnamesvar, dvnamesvar, newnc) # WRFt (T, P, PB) elif diagn == 'WRFt': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] var2 = ncobj.variables[depvars[2]][:] p0=100000. p=var1 + var2 WRFt = (var0 + 300.)*(p/p0)**(2./7.) dnamesvar = list(ncobj.variables[depvars[0]].dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'ta', WRFt, dnames, diagoutvd, newnc) # WRFtda (WRFrh, WRFt) elif diagn == 'WRFtda': ARM2 = fdef.module_definitions.arm2 ARM3 = fdef.module_definitions.arm3 gammatarh = np.log(WRFrh) + ARM2*(WRFt-273.15)/((WRFt-273.15)+ARM3) td = ARM3*gammatarh/(ARM2-gammatarh) dnamesvar = list(ncobj.variables['T'].dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'tda', td, dnames, diagoutvd, newnc) # WRFtdas (PSFC, T2, Q2) elif diagn == 'WRFtdas': ARM2 = fdef.module_definitions.arm2 ARM3 = fdef.module_definitions.arm3 var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] var2 = ncobj.variables[depvars[2]][:] dnamesvar = list(ncobj.variables[depvars[1]].dimensions) dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) rhs, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar) gammatarhs = np.log(rhs) + ARM2*(var1-273.15)/((var1-273.15)+ARM3) tdas = ARM3*gammatarhs/(ARM2-gammatarhs) + 273.15 # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'tdas', tdas, dnames, diagoutvd, newnc) # WRFua (U, V, SINALPHA, COSALPHA) to be rotated !! elif diagn == 'WRFua': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] var2 = ncobj.variables[depvars[2]][:] var3 = ncobj.variables[depvars[3]][:] # un-staggering variables if len(var0.shape) == 4: unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1] elif len(var0.shape) == 3: # Asuming sunding point (dimt, dimz, dimstgx) unstgdims = [var0.shape[0], var0.shape[1]] ua = np.zeros(tuple(unstgdims), dtype=np.float) unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float) unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float) if len(var0.shape) == 4: unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]]) unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:]) for iz in range(var0.shape[1]): ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2 dnamesvar = ['Time','bottom_top','south_north','west_east'] elif len(var0.shape) == 3: unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1]) unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1]) for iz in range(var0.shape[1]): ua[:,iz] = unstgvar0[:,iz]*var3 - unstgvar1[:,iz]*var2 dnamesvar = ['Time','bottom_top'] dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'ua', ua, dnames, diagoutvd, newnc) # WRFua (U, V, SINALPHA, COSALPHA) to be rotated !! elif diagn == 'WRFva': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] var2 = ncobj.variables[depvars[2]][:] var3 = ncobj.variables[depvars[3]][:] # un-staggering variables if len(var0.shape) == 4: unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1] elif len(var0.shape) == 3: # Asuming sunding point (dimt, dimz, dimstgx) unstgdims = [var0.shape[0], var0.shape[1]] va = np.zeros(tuple(unstgdims), dtype=np.float) unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float) unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float) if len(var0.shape) == 4: unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]]) unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:]) for iz in range(var0.shape[1]): va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3 dnamesvar = ['Time','bottom_top','south_north','west_east'] elif len(var0.shape) == 3: unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1]) unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1]) for iz in range(var0.shape[1]): va[:,iz] = unstgvar0[:,iz]*var2 + unstgvar1[:,iz]*var3 dnamesvar = ['Time','bottom_top'] dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'va', va, dnames, diagoutvd, newnc) # WRFwd (U, V, SINALPHA, COSALPHA) to be rotated !! elif diagn == 'WRFwd': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] var2 = ncobj.variables[depvars[2]][:] var3 = ncobj.variables[depvars[3]][:] # un-staggering variables if len(var0.shape) == 4: unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1] elif len(var0.shape) == 3: # Asuming sunding point (dimt, dimz, dimstgx) unstgdims = [var0.shape[0], var0.shape[1]] ua = np.zeros(tuple(unstgdims), dtype=np.float) va = np.zeros(tuple(unstgdims), dtype=np.float) unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float) unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float) if len(var0.shape) == 4: unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]]) unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:]) for iz in range(var0.shape[1]): ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2 va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3 dnamesvar = ['Time','bottom_top','south_north','west_east'] elif len(var0.shape) == 3: unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1]) unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1]) for iz in range(var0.shape[1]): ua[:,iz] = unstgvar0[:,iz]*var3 - unstgvar1[:,iz]*var2 va[:,iz] = unstgvar0[:,iz]*var2 + unstgvar1[:,iz]*var3 dnamesvar = ['Time','bottom_top'] dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) wd, dnames, dvnames = diag.compute_wd(ua, va, dnamesvar, dvnamesvar) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'wd', wd, dnames, diagoutvd, newnc) # WRFtime elif diagn == 'WRFtime': diagout = WRFtime dnamesvar = ['Time'] dvnamesvar = ['Times'] ncvar.insert_variable(ncobj, 'time', diagout, dnamesvar, dvnamesvar, newnc) # ws (U, V) elif diagn == 'ws': # un-staggering variables if len(var0.shape) == 4: unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1] elif len(var0.shape) == 3: # Asuming sunding point (dimt, dimz, dimstgx) unstgdims = [var0.shape[0], var0.shape[1]] ua = np.zeros(tuple(unstgdims), dtype=np.float) va = np.zeros(tuple(unstgdims), dtype=np.float) unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float) unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float) if len(var0.shape) == 4: unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]]) unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:]) for iz in range(var0.shape[1]): ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2 va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3 dnamesvar = ['Time','bottom_top','south_north','west_east'] elif len(var0.shape) == 3: unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1]) unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1]) for iz in range(var0.shape[1]): ua[:,iz] = unstgvar0[:,iz]*var3 - unstgvar1[:,iz]*var2 va[:,iz] = unstgvar0[:,iz]*var2 + unstgvar1[:,iz]*var3 dnamesvar = ['Time','bottom_top'] diagout = np.sqrt(unstgvar0*unstgvar0 + unstgvar1*unstgvar1) # dnamesvar = ncobj.variables[depvars[0]].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnamesvar) for nonvd in NONchkvardims: if gen.searchInlist(dvnamesvar,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'ws', diagout, dnamesvar, diagoutvd, newnc) # wss (u10, v10) elif diagn == 'wss': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] diagout = np.sqrt(var0*var0 + var1*var1) dnamesvar = ncobj.variables[depvars[0]].dimensions dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames) ncvar.insert_variable(ncobj, 'wss', diagout, dnamesvar, dvnamesvar, newnc) # WRFheight height from WRF geopotential as WRFGeop/g elif diagn == 'WRFheight': diagout = WRFgeop/grav # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'zhgt', diagout, dnames, diagoutvd, newnc) # WRFheightrel relative-height from WRF geopotential as WRFgeop(PH + PHB)/g-HGT 'WRFheightrel|PH@PHB@HGT elif diagn == 'WRFheightrel': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] var2 = ncobj.variables[depvars[2]][:] dimz = var0.shape[1] diagout = np.zeros(tuple(var0.shape), dtype=np.float) for iz in range(dimz): diagout[:,iz,:,:] = (var0[:,iz,:,:]+ var1[:,iz,:,:])/grav - var2 # Removing the nonChecking variable-dimensions from the initial list varsadd = [] diagoutvd = list(dvnames) for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'zhgtrel', diagout, dnames, diagoutvd, newnc) # WRFzmla_gen generic boundary layer hieght computation from WRF theta, QVAPOR, WRFgeop, HGT, elif diagn == 'WRFzmlagen': var0 = ncobj.variables[depvars[0]][:]+300. var1 = ncobj.variables[depvars[1]][:] dimz = var0.shape[1] var2 = WRFgeop[:,1:dimz+1,:,:]/9.8 var3 = ncobj.variables[depvars[3]][0,:,:] diagout, diagoutd, diagoutvd = diag.Forcompute_zmla_gen(var0,var1,var2,var3, \ dnames,dvnames) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'zmla', diagout, diagoutd, diagoutvd, newnc) # WRFzwind wind extrapolation at a given height using power law computation from WRF # U, V, WRFz, U10, V10, SINALPHA, COSALPHA, z=[zval] elif diagn == 'WRFzwind': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] var2 = WRFz var3 = ncobj.variables[depvars[3]][:] var4 = ncobj.variables[depvars[4]][:] var5 = ncobj.variables[depvars[5]][0,:,:] var6 = ncobj.variables[depvars[6]][0,:,:] var7 = np.float(depvars[7].split('=')[1]) # un-staggering 3D winds unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1] va = np.zeros(tuple(unstgdims), dtype=np.float) unvar0 = np.zeros(tuple(unstgdims), dtype=np.float) unvar1 = np.zeros(tuple(unstgdims), dtype=np.float) unvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]]) unvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:]) diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwind(unvar0, \ unvar1, var2, var3, var4, var5, var6, var7, dnames, dvnames) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc) ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc) # WRFzwind wind extrapolation at a given hieght using logarithmic law computation # from WRF U, V, WRFz, U10, V10, SINALPHA, COSALPHA, z=[zval] elif diagn == 'WRFzwind_log': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] var2 = WRFz var3 = ncobj.variables[depvars[3]][:] var4 = ncobj.variables[depvars[4]][:] var5 = ncobj.variables[depvars[5]][0,:,:] var6 = ncobj.variables[depvars[6]][0,:,:] var7 = np.float(depvars[7].split('=')[1]) # un-staggering 3D winds unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1] va = np.zeros(tuple(unstgdims), dtype=np.float) unvar0 = np.zeros(tuple(unstgdims), dtype=np.float) unvar1 = np.zeros(tuple(unstgdims), dtype=np.float) unvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]]) unvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:]) diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwind_log(unvar0, \ unvar1, var2, var3, var4, var5, var6, var7, dnames, dvnames) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc) ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc) # WRFzwindMO wind extrapolation at a given height computation using Monin-Obukhow # theory from WRF UST, ZNT, RMOL, U10, V10, SINALPHA, COSALPHA, z=[zval] # NOTE: only useful for [zval] < 80. m elif diagn == 'WRFzwindMO': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] var2 = ncobj.variables[depvars[2]][:] var3 = ncobj.variables[depvars[3]][:] var4 = ncobj.variables[depvars[4]][:] var5 = ncobj.variables[depvars[5]][0,:,:] var6 = ncobj.variables[depvars[6]][0,:,:] var7 = np.float(depvars[7].split('=')[1]) diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwindMO(var0, var1,\ var2, var3, var4, var5, var6, var7, dnames, dvnames) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc) ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc) # zmla_gen generic boundary layer hieght computation from generic 2D-file theta, qv, zg, orog, elif diagn == 'zmlagen': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] dimz = var0.shape[1] var2 = ncobj.variables[depvars[2]][:] var3 = ncobj.variables[depvars[3]][:] diagout, diagoutd, diagoutvd = diag.Forcompute_zmla_gen(var0,var1,var2,var3, \ dnames,dvnames) # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'zmla', diagout, diagoutd, diagoutvd, newnc) # zmla_genUWsnd generic boundary layer hieght computation from UWyoming sounding file theta, qv, zg elif diagn == 'zmlagenUWsnd': var0 = ncobj.variables[depvars[0]][:] var1 = ncobj.variables[depvars[1]][:] dimz = var0.shape[1] var2 = ncobj.variables[depvars[2]][:]/9.8 var3 = ncobj.getncattr('Station_elevation') diagout, diagoutd, diagoutvd = diag.Forcompute_zmla_gen(var0,var1,var2,var3, \ dnames,dvnames) # No need to remove topography height diagout = diagout + var3 # Removing the nonChecking variable-dimensions from the initial list varsadd = [] for nonvd in NONchkvardims: if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd) varsadd.append(nonvd) ncvar.insert_variable(ncobj, 'zmla', diagout, diagoutd, diagoutvd, newnc) else: print errormsg print ' ' + main + ": diagnostic '" + diagn + "' not ready!!!" print ' available diagnostics: ', availdiags quit(-1) newnc.sync() # Adding that additional variables required to compute some diagnostics which # where not in the original file print ' adding additional variables...' for vadd in varsadd: if not gen.searchInlist(newnc.variables.keys(),vadd) and \ dictcompvars.has_key(vadd): attrs = dictcompvars[vadd] vvn = attrs['name'] if not gen.searchInlist(newnc.variables.keys(), vvn): iidvn = dvnames.index(vadd) dnn = dnames[iidvn] if vadd == 'WRFtime': dvarvals = WRFtime[:] newvar = newnc.createVariable(vvn, 'f8', (dnn)) newvar[:] = dvarvals for attn in attrs.keys(): if attn != 'name': attv = attrs[attn] ncvar.set_attribute(newvar, attn, attv) # end of diagnostics # Global attributes ## ncvar.add_global_PyNCplot(newnc, main, None, '2.0') gorigattrs = ncobj.ncattrs() for attr in gorigattrs: attrv = ncobj.getncattr(attr) atvar = ncvar.set_attribute(newnc, attr, attrv) ncobj.close() newnc.close() print '\n' + main + ': successfull writting of diagnostics file "' + ofile + '" !!!'