Context Navigation

source: lmdz_wrf/trunk/tools/diagnostics.py @ 2659

Last change on this file since 2659 was 2655, checked in by lfita, 6 years ago

Adding:

`front_R04': generalization of front calculation from Rodrigues, 2004

File size: 92.1 KB

Rev	Line
[1675]	1	# Python script to comput diagnostics
[1908]	2	# From L. Fita work in different places: CCRC (Australia), LMD (France)
	3	# More information at: http://www.xn--llusfb-5va.cat/python/PyNCplot
	4	#
	5	# pyNCplot and its component nc_var.py comes with ABSOLUTELY NO WARRANTY.
	6	# This work is licendes under a Creative Commons
	7	# Attribution-ShareAlike 4.0 International License (http://creativecommons.org/licenses/by-sa/4.0)
	8	#
	9	# L. Fita, CIMA. CONICET-UBA, CNRS UMI-IFAECI, C.A. Buenos Aires, Argentina
[365]	10	# File diagnostics.inf provides the combination of variables to get the desired diagnostic
[772]	11	# To be used with module_ForDiagnostics.F90, module_ForDiagnosticsVars.F90, module_generic.F90
[1150]	12	# foudre: f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90 >& run_f2py.log
	13	# 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
[1149]	14
[2095]	15	## 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
[413]	16	## e.g. # diagnostics.py -f /home/lluis/PY/diagnostics.inf -d variable_combo -v WRFprc
[365]	17
[1675]	18	# Available general pupose diagnostics (model independent) providing (varv1, varv2, ..., dimns, dimvns)
	19	# compute_accum: Function to compute the accumulation of a variable
	20	# compute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following
	21	# newmicro.F90 from LMDZ compute_clt(cldfra, pres, dimns, dimvns)
	22	# compute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ
	23	# compute_clivi: Function to compute cloud-ice water path (clivi)
	24	# compute_clwvl: Function to compute condensed water path (clwvl)
	25	# compute_deaccum: Function to compute the deaccumulation of a variable
	26	# compute_mslp: Function to compute mslp: mean sea level pressure following p_interp.F90 from WRF
	27	# compute_OMEGAw: Function to transform OMEGA [Pas-1] to velocities [ms-1]
	28	# compute_prw: Function to compute water vapour path (prw)
[2209]	29	# compute_range_faces: Function to compute faces [uphill, valley, downhill] of sections of a mountain
	30	# range, along a given face
[1675]	31	# compute_rh: Function to compute relative humidity following 'Tetens' equation (T,P) ...'
	32	# compute_td: Function to compute the dew point temperature
	33	# compute_turbulence: Function to compute the rubulence term of the Taylor's decomposition ...'
	34	# compute_wds: Function to compute the wind direction
	35	# compute_wss: Function to compute the wind speed
	36	# compute_WRFuava: Function to compute geographical rotated WRF 3D winds
	37	# compute_WRFuasvas: Fucntion to compute geographical rotated WRF 2-meter winds
	38	# derivate_centered: Function to compute the centered derivate of a given field
	39	# def Forcompute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ via Fortran subroutine
	40	# Forcompute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ via a Fortran module
[1758]	41	# Forcompute_psl_ptarget: Function to compute the sea-level pressure following target_pressure value found in `p_interp.F'
[1675]	42
	43	# Others just providing variable values
	44	# var_cllmh: Fcuntion to compute cllmh on a 1D column
	45	# var_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ using 1D vertical column values
	46	# var_mslp: Fcuntion to compute mean sea-level pressure
	47	# var_virtualTemp: This function returns virtual temperature in K,
	48	# var_WRFtime: Function to copmute CFtimes from WRFtime variable
	49	# rotational_z: z-component of the rotatinoal of horizontal vectorial field
	50	# turbulence_var: Function to compute the Taylor's decomposition turbulence term from a a given variable
[2138]	51	# timeoverthres: When a given variable [varname] overpass a given [value]. Being [CFvarn] the name of the diagnostics in
	52	# variables_values.dat
[2141]	53	# timemax ([varname], time). When a given variable [varname] got its maximum
[1675]	54
[365]	55	from optparse import OptionParser
	56	import numpy as np
	57	from netCDF4 import Dataset as NetCDFFile
	58	import os
	59	import re
	60	import nc_var_tools as ncvar
[756]	61	import generic_tools as gen
[654]	62	import datetime as dtime
[1163]	63	import module_ForDiag as fdin
[1942]	64	import module_ForDef as fdef
[1675]	65	import diag_tools as diag
[365]	66
	67	main = 'diagnostics.py'
	68	errormsg = 'ERROR -- error -- ERROR -- error'
	69	warnmsg = 'WARNING -- warning -- WARNING -- warning'
	70
[654]	71	# Constants
	72	grav = 9.81
	73
[365]	74
	75	####### ###### ##### #### ### ## #
	76	comboinf="\nIF -d 'variable_combo', provides information of the combination to obtain -v [varn] with the ASCII file with the combinations as -f [combofile]"
	77
	78	parser = OptionParser()
	79	parser.add_option("-f", "--netCDF_file", dest="ncfile", help="file to use", metavar="FILE")
	80	parser.add_option("-d", "--dimensions", dest="dimns",
[1761]	81	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,
[365]	82	metavar="LABELS")
	83	parser.add_option("-v", "--variables", dest="varns",
	84	help=" [varn1]\|[var11]@[...[varN1]],[...,[varnM]\|[var1M]@[...[varLM]]] ',' list of variables to compute [varnK] and its necessary ones [var1K]...[varPK]", metavar="VALUES")
	85
	86	(opts, args) = parser.parse_args()
	87
	88	####### #######
	89	## MAIN
	90	#######
[2100]	91	availdiags = ['ACRAINTOT', 'accum', 'clt', 'cllmh', 'convini', 'deaccum', 'fog_K84', \
[2643]	92	'fog_RUC', 'front_R04', 'rhs_tas_tds', 'LMDZrh', 'mslp', 'OMEGAw', 'RAINTOT', \
	93	'range_faces', \
[2390]	94	'rvors', 'td', 'timemax', 'timeoverthres', 'turbulence', 'tws', 'uavaFROMwswd', \
[2274]	95	'WRFbnds', 'WRFcape_afwa', 'WRFclivi', 'WRFclwvi', 'WRF_denszint', 'WRFgeop', \
[2206]	96	'WRFmrso', 'WRFmrsos', 'WRFpotevap_orPM', 'WRFp', 'WRFpsl_ecmwf', \
[1762]	97	'WRFpsl_ptarget', 'WRFrvors', 'WRFslw', 'ws', 'wds', 'wss', 'WRFheight', \
[2390]	98	'WRFheightrel', 'WRFtda', 'WRFtdas', 'WRFtws', 'WRFua', 'WRFva', 'WRFzwind', \
[2621]	99	'WRFzwind_log', 'WRFzwindMO', 'zmlagen', 'zmlagenUWsnd']
[365]	100
[649]	101	methods = ['accum', 'deaccum']
	102
[365]	103	# Variables not to check
[2277]	104	NONcheckingvars = ['accum', 'cllmh', 'deaccum', 'face', 'LONLATdxdy', \
[2655]	105	'params', 'reglonlatbnds', 'TSrhs', 'TStd', 'TSwds', 'TSwss', \
[2274]	106	'WRFbils', 'WRFbnds', \
[2257]	107	'WRFclivi', 'WRFclwvi', 'WRFdens', 'WRFdx', 'WRFdxdy', 'WRFdxdywps', 'WRFdy', \
	108	'WRFgeop', 'WRFp', 'WRFtd', \
[1809]	109	'WRFpos', 'WRFprc', 'WRFprls', 'WRFrh', 'LMDZrh', 'LMDZrhs', \
[1806]	110	'WRFrhs', 'WRFrvors', \
[2215]	111	'WRFt', 'WRFtime', 'WRFua', 'WRFva', 'WRFwds', 'WRFwss', 'WRFheight', 'WRFz', \
	112	'WRFzg']
[365]	113
[1809]	114	# diagnostics not to check their dependeny
[2138]	115	NONcheckdepvars = ['accum', 'deaccum', 'timeoverthres', 'WRF_denszint', \
	116	'WRFzwind_log', 'WRFzwind', 'WRFzwindMO']
[1809]	117
[1351]	118	NONchkvardims = ['WRFtime']
	119
[365]	120	ofile = 'diagnostics.nc'
	121
	122	dimns = opts.dimns
	123	varns = opts.varns
	124
	125	# Special method. knowing variable combination
	126	##
	127	if opts.dimns == 'variable_combo':
	128	print warnmsg
	129	print ' ' + main + ': knowing variable combination !!!'
	130	combination = variable_combo(opts.varns,opts.ncfile)
	131	print ' COMBO: ' + combination
	132	quit(-1)
	133
[1833]	134	if opts.ncfile is None:
	135	print errormsg
	136	print ' ' + main + ": No file provided !!"
	137	print ' is mandatory to provide a file -f [filename]'
	138	quit(-1)
	139
	140	if opts.dimns is None:
	141	print errormsg
	142	print ' ' + main + ": No description of dimensions are provided !!"
	143	print ' is mandatory to provide description of dimensions as -d [dimn]@[vardimname],... '
	144	quit(-1)
	145
	146	if opts.varns is None:
	147	print errormsg
	148	print ' ' + main + ": No variable to diagnose is provided !!"
	149	print ' is mandatory to provide a variable to diagnose as -v [diagn]\|[varn1]@... '
	150	quit(-1)
	151
[365]	152	if not os.path.isfile(opts.ncfile):
	153	print errormsg
	154	print ' ' + main + ": file '" + opts.ncfile + "' does not exist !!"
	155	quit(-1)
	156
	157	ncobj = NetCDFFile(opts.ncfile, 'r')
	158
[1351]	159	# Looking for specific variables that might be use in more than one diagnostic
	160	WRFgeop_compute = False
	161	WRFp_compute = False
	162	WRFt_compute = False
	163	WRFrh_compute = False
	164	WRFght_compute = False
	165	WRFdens_compute = False
	166	WRFpos_compute = False
	167	WRFtime_compute = False
[1777]	168	WRFz_compute = False
[2655]	169	WRFdx_compute = False
	170	WRFdx_compute = False
[2215]	171	WRFdxdy_compute = False
[2257]	172	WRFdxdywps_compute = False
[2222]	173	LONLATdxdy_compute = False
[1351]	174
[365]	175	# File creation
	176	newnc = NetCDFFile(ofile,'w')
	177
	178	# dimensions
	179	dimvalues = dimns.split(',')
	180	dnames = []
	181	dvnames = []
	182
	183	for dimval in dimvalues:
[1351]	184	dn = dimval.split('@')[0]
	185	dnv = dimval.split('@')[1]
	186	dnames.append(dn)
	187	dvnames.append(dnv)
	188	# Is there any dimension-variable which should be computed?
	189	if dnv == 'WRFgeop':WRFgeop_compute = True
	190	if dnv == 'WRFp': WRFp_compute = True
	191	if dnv == 'WRFt': WRFt_compute = True
	192	if dnv == 'WRFrh': WRFrh_compute = True
	193	if dnv == 'WRFght': WRFght_compute = True
	194	if dnv == 'WRFdens': WRFdens_compute = True
	195	if dnv == 'WRFpos': WRFpos_compute = True
	196	if dnv == 'WRFtime': WRFtime_compute = True
[1777]	197	if dnv == 'WRFz':WRFz_compute = True
[2655]	198	if dnv == 'WRFdx':WRFdx_compute = True
	199	if dnv == 'WRFdy':WRFdy_compute = True
[2215]	200	if dnv == 'WRFdxdy':WRFdxdy_compute = True
[2257]	201	if dnv == 'WRFdxdywps':WRFdxdywps_compute = True
[2222]	202	if dnv == 'LONLATdxdy':LONLATdxdy_compute = True
[365]	203
	204	# diagnostics to compute
	205	diags = varns.split(',')
	206	Ndiags = len(diags)
	207
	208	for idiag in range(Ndiags):
	209	if diags[idiag].split('\|')[1].find('@') == -1:
	210	depvars = diags[idiag].split('\|')[1]
[654]	211	if depvars == 'WRFgeop':WRFgeop_compute = True
[365]	212	if depvars == 'WRFp': WRFp_compute = True
	213	if depvars == 'WRFt': WRFt_compute = True
	214	if depvars == 'WRFrh': WRFrh_compute = True
	215	if depvars == 'WRFght': WRFght_compute = True
	216	if depvars == 'WRFdens': WRFdens_compute = True
	217	if depvars == 'WRFpos': WRFpos_compute = True
[654]	218	if depvars == 'WRFtime': WRFtime_compute = True
[1777]	219	if depvars == 'WRFz': WRFz_compute = True
[365]	220	else:
	221	depvars = diags[idiag].split('\|')[1].split('@')
[756]	222	if gen.searchInlist(depvars, 'WRFgeop'): WRFgeop_compute = True
	223	if gen.searchInlist(depvars, 'WRFp'): WRFp_compute = True
	224	if gen.searchInlist(depvars, 'WRFt'): WRFt_compute = True
	225	if gen.searchInlist(depvars, 'WRFrh'): WRFrh_compute = True
	226	if gen.searchInlist(depvars, 'WRFght'): WRFght_compute = True
	227	if gen.searchInlist(depvars, 'WRFdens'): WRFdens_compute = True
	228	if gen.searchInlist(depvars, 'WRFpos'): WRFpos_compute = True
	229	if gen.searchInlist(depvars, 'WRFtime'): WRFtime_compute = True
[1777]	230	if gen.searchInlist(depvars, 'WRFz'): WRFz_compute = True
[2655]	231	if gen.searchInlist(depvars, 'WRFdx'): WRFdx_compute = True
	232	if gen.searchInlist(depvars, 'WRFdy'): WRFdy_compute = True
[2215]	233	if gen.searchInlist(depvars, 'WRFdxdy'): WRFdxdy_compute = True
[2257]	234	if gen.searchInlist(depvars, 'WRFdxdywps'): WRFdxdywps_compute = True
[2222]	235	if gen.searchInlist(depvars, 'LONLATdxdy'): LONLATdxdy_compute = True
[365]	236
[1351]	237	# Dictionary with the new computed variables to be able to add them
	238	dictcompvars = {}
[654]	239	if WRFgeop_compute:
	240	print ' ' + main + ': Retrieving geopotential value from WRF as PH + PHB'
	241	dimv = ncobj.variables['PH'].shape
	242	WRFgeop = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
	243
[1351]	244	# Attributes of the variable
[1412]	245	Vvals = gen.variables_values('WRFgeop')
[1351]	246	dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
	247	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
	248
[365]	249	if WRFp_compute:
	250	print ' ' + main + ': Retrieving pressure value from WRF as P + PB'
	251	dimv = ncobj.variables['P'].shape
	252	WRFp = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	253
[1351]	254	# Attributes of the variable
	255	Vvals = gen.variables_values('WRFp')
	256	dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
	257	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
	258
[365]	259	if WRFght_compute:
	260	print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...'
	261	WRFght = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
	262
[1351]	263	# Attributes of the variable
	264	Vvals = gen.variables_values('WRFght')
	265	dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
	266	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
	267
[365]	268	if WRFrh_compute:
	269	print ' ' + main + ": computing relative humidity from WRF as 'Tetens'" + \
	270	' equation (T,P) ...'
	271	p0=100000.
	272	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	273	tk = (ncobj.variables['T'][:] + 300.)(p/p0)*(2./7.)
	274	qv = ncobj.variables['QVAPOR'][:]
	275
	276	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
	277	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
	278
	279	WRFrh = qv/data2
	280
[1351]	281	# Attributes of the variable
	282	Vvals = gen.variables_values('WRFrh')
	283	dictcompvars['WRFrh'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
	284	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
	285
[365]	286	if WRFt_compute:
	287	print ' ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
	288	p0=100000.
	289	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	290
	291	WRFt = (ncobj.variables['T'][:] + 300.)(p/p0)*(2./7.)
	292
[1351]	293	# Attributes of the variable
	294	Vvals = gen.variables_values('WRFt')
	295	dictcompvars['WRFt'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
	296	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
	297
[365]	298	if WRFdens_compute:
	299	print ' ' + main + ': computing air density from WRF as ((MU + MUB) * ' + \
	300	'DNW)/g ...'
	301
	302	# Just we need in in absolute values: Size of the central grid cell
	303	## dxval = ncobj.getncattr('DX')
	304	## dyval = ncobj.getncattr('DY')
	305	## mapfac = ncobj.variables['MAPFAC_M'][:]
	306	## area = dxvaldyvalmapfac
	307
	308	mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
	309	dnw = ncobj.variables['DNW'][:]
	310
	311	WRFdens = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]), \
	312	dtype=np.float)
	313	levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)
	314
	315	for it in range(mu.shape[0]):
	316	for iz in range(dnw.shape[1]):
	317	levval.fill(np.abs(dnw[it,iz]))
	318	WRFdens[it,iz,:,:] = levval
	319	WRFdens[it,iz,:,:] = mu[it,:,:]*WRFdens[it,iz,:,:]/grav
	320
[1351]	321	# Attributes of the variable
	322	Vvals = gen.variables_values('WRFdens')
	323	dictcompvars['WRFdens'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
	324	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
	325
[365]	326	if WRFpos_compute:
	327	# WRF positions from the lowest-leftest corner of the matrix
	328	print ' ' + main + ': computing position from MAPFAC_M as sqrt(DYj*2 + ' + \
	329	'DXx2)MAPFAC_M ...'
	330
	331	mapfac = ncobj.variables['MAPFAC_M'][:]
	332
	333	distx = np.float(ncobj.getncattr('DX'))
	334	disty = np.float(ncobj.getncattr('DY'))
	335
	336	print 'distx:',distx,'disty:',disty
	337
	338	dx = mapfac.shape[2]
	339	dy = mapfac.shape[1]
	340	dt = mapfac.shape[0]
	341
	342	WRFpos = np.zeros((dt, dy, dx), dtype=np.float)
	343
	344	for i in range(1,dx):
	345	WRFpos[0,0,i] = distx*i/mapfac[0,0,i]
	346	for j in range(1,dy):
	347	i=0
	348	WRFpos[0,j,i] = WRFpos[0,j-1,i] + disty/mapfac[0,j,i]
	349	for i in range(1,dx):
	350	# WRFpos[0,j,i] = np.sqrt((distyj)2. + (distxi)**2.)/mapfac[0,j,i]
	351	# WRFpos[0,j,i] = np.sqrt((distyj)2. + (distxi)**2.)
	352	WRFpos[0,j,i] = WRFpos[0,j,i-1] + distx/mapfac[0,j,i]
	353
	354	for it in range(1,dt):
	355	WRFpos[it,:,:] = WRFpos[0,:,:]
	356
[654]	357	if WRFtime_compute:
	358	print ' ' + main + ': computing time from WRF as CFtime(Times) ...'
	359
	360	refdate='19491201000000'
	361	tunitsval='minutes'
	362
	363	timeobj = ncobj.variables['Times']
	364	timewrfv = timeobj[:]
	365
	366	yrref=refdate[0:4]
	367	monref=refdate[4:6]
	368	dayref=refdate[6:8]
	369	horref=refdate[8:10]
	370	minref=refdate[10:12]
	371	secref=refdate[12:14]
	372
	373	refdateS = yrref + '-' + monref + '-' + dayref + ' ' + horref + ':' + minref + \
	374	':' + secref
	375
[2065]	376
	377	if len(timeobj.shape) == 2:
	378	dt = timeobj.shape[0]
	379	else:
	380	dt = 1
[654]	381	WRFtime = np.zeros((dt), dtype=np.float)
	382
[2065]	383	if len(timeobj.shape) == 2:
	384	for it in range(dt):
	385	wrfdates = gen.datetimeStr_conversion(timewrfv[it,:],'WRFdatetime', \
	386	'matYmdHMS')
	387	WRFtime[it] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
	388	else:
	389	wrfdates = gen.datetimeStr_conversion(timewrfv[:],'WRFdatetime', \
	390	'matYmdHMS')
[2072]	391	WRFtime[0] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
[654]	392
	393	tunits = tunitsval + ' since ' + refdateS
	394
[1351]	395	# Attributes of the variable
	396	dictcompvars['WRFtime'] = {'name': 'time', 'standard_name': 'time', \
	397	'long_name': 'time', 'units': tunits, 'calendar': 'gregorian'}
	398
[1777]	399	if WRFz_compute:
	400	print ' ' + main + ': Retrieving z: height above surface value from WRF as ' + \
	401	'unstagger(PH + PHB)/9.8-hgt'
	402	dimv = ncobj.variables['PH'].shape
	403	WRFzg = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/9.8
	404
	405	unzgd = (dimv[0], dimv[1]-1, dimv[2], dimv[3])
	406	unzg = np.zeros(unzgd, dtype=np.float)
	407	unzg = 0.5*(WRFzg[:,0:dimv[1]-1,:,:] + WRFzg[:,1:dimv[1],:,:])
	408
	409	WRFz = np.zeros(unzgd, dtype=np.float)
	410	for iz in range(dimv[1]-1):
	411	WRFz[:,iz,:,:] = unzg[:,iz,:,:] - ncobj.variables['HGT'][:]
	412
	413	# Attributes of the variable
	414	Vvals = gen.variables_values('WRFz')
[2215]	415	dictcompvars['WRFz'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
[1777]	416	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
	417
[2655]	418	if WRFdxdy_compute or WRFdx_compute or WRFdy_copute:
[2215]	419	print ' ' + main + ': Retrieving dxdy: real distance between grid points ' + \
	420	'from WRF as dx=(XLONG(i+1)-XLONG(i))DX/MAPFAC_M, dy=(XLAT(j+1)-XLAT(i))DY/'+\
	421	'MAPFAC_M, ds=sqrt(dx2+dy2)'
	422	dimv = ncobj.variables['XLONG'].shape
	423	WRFlon = ncobj.variables['XLONG'][0,:,:]
	424	WRFlat = ncobj.variables['XLAT'][0,:,:]
	425	WRFmapfac_m = ncobj.variables['MAPFAC_M'][0,:,:]
	426	DX = ncobj.DX
	427	DY = ncobj.DY
	428
	429	dimx = dimv[2]
	430	dimy = dimv[1]
	431
	432	WRFdx = np.zeros((dimy,dimx), dtype=np.float)
	433	WRFdy = np.zeros((dimy,dimx), dtype=np.float)
	434
	435	WRFdx[:,0:dimx-1]=(WRFlon[:,1:dimx]-WRFlon[:,0:dimx-1])*DX/WRFmapfac_m[:,0:dimx-1]
	436	WRFdy[0:dimy-1,:]=(WRFlat[1:dimy,:]-WRFlat[0:dimy-1,:])*DY/WRFmapfac_m[0:dimy-1,:]
	437	WRFds = np.sqrt(WRFdx2 + WRFdy2)
	438
	439	# Attributes of the variable
	440	Vvals = gen.variables_values('WRFdx')
	441	dictcompvars['WRFdx'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
	442	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
	443	Vvals = gen.variables_values('WRFdy')
	444	dictcompvars['WRFdy'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
	445	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
	446	Vvals = gen.variables_values('WRFds')
[2222]	447	dictcompvars['WRFds'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
	448	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
	449
[2257]	450	if WRFdxdywps_compute:
	451	print ' ' + main + ': Retrieving dxdy: real distance between grid points ' + \
	452	'from wpsWRF as dx=(XLONG_M(i+1)-XLONG_M(i))*DX/MAPFAC_M, ' + \
	453	'dy=(XLAT_M(j+1)-XLAT_M(i))DY/MAPFAC_M, ds=sqrt(dx2+dy*2)'
	454	dimv = ncobj.variables['XLONG_M'].shape
	455	WRFlon = ncobj.variables['XLONG_M'][0,:,:]
	456	WRFlat = ncobj.variables['XLAT_M'][0,:,:]
	457	WRFmapfac_m = ncobj.variables['MAPFAC_M'][0,:,:]
	458	DX = ncobj.DX
	459	DY = ncobj.DY
	460
	461	dimx = dimv[2]
	462	dimy = dimv[1]
	463
	464	WRFdx = np.zeros((dimy,dimx), dtype=np.float)
	465	WRFdy = np.zeros((dimy,dimx), dtype=np.float)
	466
	467	WRFdx[:,0:dimx-1]=(WRFlon[:,1:dimx]-WRFlon[:,0:dimx-1])*DX/WRFmapfac_m[:,0:dimx-1]
	468	WRFdy[0:dimy-1,:]=(WRFlat[1:dimy,:]-WRFlat[0:dimy-1,:])*DY/WRFmapfac_m[0:dimy-1,:]
	469	WRFds = np.sqrt(WRFdx2 + WRFdy2)
	470
	471	# Attributes of the variable
	472	Vvals = gen.variables_values('WRFdx')
	473	dictcompvars['WRFdx'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
	474	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
	475	Vvals = gen.variables_values('WRFdy')
	476	dictcompvars['WRFdy'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
	477	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
	478	Vvals = gen.variables_values('WRFds')
	479	dictcompvars['WRFds'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
	480	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
	481
[2222]	482	if LONLATdxdy_compute:
	483	print ' ' + main + ': Retrieving dxdy: real distance between grid points ' + \
	484	'from a regular lonlat projection as dx=(lon[i+1]-lon[i])raddegRearth*' + \
	485	'cos(abs(lat[i])); dy=(lat[j+1]-lat[i])raddegRearth; ds=sqrt(dx2+dy2); '+\
	486	'raddeg = pi/180; Rearth=6370.0e03'
	487	dimv = ncobj.variables['lon'].shape
	488	lon = ncobj.variables['lon'][:]
	489	lat = ncobj.variables['lat'][:]
	490
	491	WRFlon, WRFlat = gen.lonlat2D(lon,lat)
	492
	493	dimx = WRFlon.shape[1]
	494	dimy = WRFlon.shape[0]
	495
	496	WRFdx = np.zeros((dimy,dimx), dtype=np.float)
	497	WRFdy = np.zeros((dimy,dimx), dtype=np.float)
	498
	499	raddeg = np.pi/180.
	500
	501	Rearth = fdef.module_definitions.earthradii
	502
	503	WRFdx[:,0:dimx-1]=(WRFlon[:,1:dimx]-WRFlon[:,0:dimx-1])raddegRearth* \
	504	np.cos(np.abs(WRFlat[:,0:dimx-1]*raddeg))
	505	WRFdy[0:dimy-1,:]=(WRFlat[1:dimy,:]-WRFlat[0:dimy-1,:])raddegRearth
	506	WRFds = np.sqrt(WRFdx2 + WRFdy2)
	507
	508	# Attributes of the variable
	509	Vvals = gen.variables_values('WRFdx')
	510	dictcompvars['WRFdx'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
	511	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
	512	Vvals = gen.variables_values('WRFdy')
	513	dictcompvars['WRFdy'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
	514	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
	515	Vvals = gen.variables_values('WRFds')
[2215]	516	dictcompvars['WRFds'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
	517	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
	518
[365]	519	### ## #
	520	# Going for the diagnostics
	521	### ## #
	522	print ' ' + main + ' ...'
[1404]	523	varsadd = []
[365]	524
[2390]	525	Varns = ncobj.variables.keys()
	526	Varns.sort()
	527
[365]	528	for idiag in range(Ndiags):
	529	print ' diagnostic:',diags[idiag]
[1758]	530	diagn = diags[idiag].split('\|')[0]
[365]	531	depvars = diags[idiag].split('\|')[1].split('@')
[1809]	532	if not gen.searchInlist(NONcheckdepvars, diagn):
	533	if diags[idiag].split('\|')[1].find('@') != -1:
	534	depvars = diags[idiag].split('\|')[1].split('@')
	535	if depvars[0] == 'deaccum': diagn='deaccum'
	536	if depvars[0] == 'accum': diagn='accum'
	537	for depv in depvars:
[2212]	538	# Checking without extra arguments of a depending variable (':', separated)
	539	if depv.find(':') != -1: depv=depv.split(':')[0]
[1809]	540	if not ncobj.variables.has_key(depv) and not \
	541	gen.searchInlist(NONcheckingvars, depv) and \
	542	not gen.searchInlist(methods, depv) and not depvars[0] == 'deaccum'\
	543	and not depvars[0] == 'accum' and not depv[0:2] == 'z=':
	544	print errormsg
	545	print ' ' + main + ": file '" + opts.ncfile + \
	546	"' does not have variable '" + depv + "' !!"
[2390]	547	print ' available ones:', Varns
[1809]	548	quit(-1)
	549	else:
	550	depvars = diags[idiag].split('\|')[1]
	551	if not ncobj.variables.has_key(depvars) and not \
	552	gen.searchInlist(NONcheckingvars, depvars) and \
	553	not gen.searchInlist(methods, depvars):
[365]	554	print errormsg
	555	print ' ' + main + ": file '" + opts.ncfile + \
[1809]	556	"' does not have variable '" + depvars + "' !!"
[2390]	557	print ' available ones:', Varns
[365]	558	quit(-1)
	559
[1758]	560	print "\n Computing '" + diagn + "' from: ", depvars, '...'
[365]	561
[2095]	562	# acraintot: accumulated total precipitation from WRF RAINC, RAINNC, RAINSH
[1758]	563	if diagn == 'ACRAINTOT':
[365]	564
	565	var0 = ncobj.variables[depvars[0]]
	566	var1 = ncobj.variables[depvars[1]]
[2095]	567	var2 = ncobj.variables[depvars[2]]
[365]	568
[2095]	569	diagout = var0[:] + var1[:] + var2[:]
	570
[365]	571	dnamesvar = var0.dimensions
	572	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	573
[1647]	574	# Removing the nonChecking variable-dimensions from the initial list
	575	varsadd = []
	576	for nonvd in NONchkvardims:
	577	if gen.searchInlist(dvnamesvar,nonvd): dvnamesvar.remove(nonvd)
	578	varsadd.append(nonvd)
	579
[649]	580	ncvar.insert_variable(ncobj, 'pracc', diagout, dnamesvar, dvnamesvar, newnc)
[365]	581
[649]	582	# accum: acumulation of any variable as (Variable, time [as [tunits]
	583	# from/since ....], newvarname)
[1758]	584	elif diagn == 'accum':
[649]	585
	586	var0 = ncobj.variables[depvars[0]]
	587	var1 = ncobj.variables[depvars[1]]
	588
	589	dnamesvar = var0.dimensions
	590	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	591
[1675]	592	diagout, diagoutd, diagoutvd = diag.compute_accum(var0,dnamesvar,dvnamesvar)
[1825]	593	# Removing the nonChecking variable-dimensions from the initial list
	594	varsadd = []
	595	diagoutvd = list(dvnames)
	596	for nonvd in NONchkvardims:
	597	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	598	varsadd.append(nonvd)
[649]	599
	600	CFvarn = ncvar.variables_values(depvars[0])[0]
	601
	602	# Removing the flux
	603	if depvars[1] == 'XTIME':
	604	dtimeunits = var1.getncattr('description')
	605	tunits = dtimeunits.split(' ')[0]
	606	else:
	607	dtimeunits = var1.getncattr('units')
	608	tunits = dtimeunits.split(' ')[0]
	609
[1825]	610	dtime = (var1[1] - var1[0])*diag.timeunits_seconds(tunits)
[649]	611
	612	ncvar.insert_variable(ncobj, CFvarn + 'acc', diagout*dtime, diagoutd, diagoutvd, newnc)
	613
[365]	614	# cllmh with cldfra, pres
[1758]	615	elif diagn == 'cllmh':
[365]	616
	617	var0 = ncobj.variables[depvars[0]]
	618	if depvars[1] == 'WRFp':
	619	var1 = WRFp
	620	else:
	621	var01 = ncobj.variables[depvars[1]]
	622	if len(size(var1.shape)) < len(size(var0.shape)):
	623	var1 = np.brodcast_arrays(var01,var0)[0]
	624	else:
	625	var1 = var01
	626
[1675]	627	diagout, diagoutd, diagoutvd = diag.Forcompute_cllmh(var0,var1,dnames,dvnames)
[772]	628
[1351]	629	# Removing the nonChecking variable-dimensions from the initial list
	630	varsadd = []
	631	for nonvd in NONchkvardims:
	632	if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
	633	varsadd.append(nonvd)
	634
[365]	635	ncvar.insert_variable(ncobj, 'cll', diagout[0,:], diagoutd, diagoutvd, newnc)
	636	ncvar.insert_variable(ncobj, 'clm', diagout[1,:], diagoutd, diagoutvd, newnc)
	637	ncvar.insert_variable(ncobj, 'clh', diagout[2,:], diagoutd, diagoutvd, newnc)
	638
	639	# clt with cldfra
[1758]	640	elif diagn == 'clt':
[365]	641
	642	var0 = ncobj.variables[depvars]
[1675]	643	diagout, diagoutd, diagoutvd = diag.Forcompute_clt(var0,dnames,dvnames)
[1351]	644
	645	# Removing the nonChecking variable-dimensions from the initial list
	646	varsadd = []
	647	for nonvd in NONchkvardims:
	648	if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
	649	varsadd.append(nonvd)
	650
[365]	651	ncvar.insert_variable(ncobj, 'clt', diagout, diagoutd, diagoutvd, newnc)
	652
[2100]	653	# convini (pr, time)
	654	elif diagn == 'convini':
	655
	656	var0 = ncobj.variables[depvars[0]][:]
	657	var1 = ncobj.variables[depvars[1]][:]
	658	otime = ncobj.variables[depvars[1]]
	659
	660	dnamesvar = ncobj.variables[depvars[0]].dimensions
	661	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	662
	663	diagout, diagoutd, diagoutvd = diag.var_convini(var0, var1, dnames, dvnames)
	664
	665	ncvar.insert_variable(ncobj, 'convini', diagout, diagoutd, diagoutvd, newnc, \
	666	fill=gen.fillValueF)
	667	# Getting the right units
	668	ovar = newnc.variables['convini']
	669	if gen.searchInlist(otime.ncattrs(), 'units'):
	670	tunits = otime.getncattr('units')
	671	ncvar.set_attribute(ovar, 'units', tunits)
	672	newnc.sync()
	673
[365]	674	# deaccum: deacumulation of any variable as (Variable, time [as [tunits]
	675	# from/since ....], newvarname)
[1758]	676	elif diagn == 'deaccum':
[365]	677
[1825]	678	var0 = ncobj.variables[depvars[0]]
	679	var1 = ncobj.variables[depvars[1]]
[365]	680
	681	dnamesvar = var0.dimensions
	682	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	683
[1675]	684	diagout, diagoutd, diagoutvd = diag.compute_deaccum(var0,dnamesvar,dvnamesvar)
[1825]	685	# Removing the nonChecking variable-dimensions from the initial list
	686	varsadd = []
	687	diagoutvd = list(dvnames)
	688	for nonvd in NONchkvardims:
	689	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	690	varsadd.append(nonvd)
[365]	691
	692	# Transforming to a flux
[1825]	693	if depvars[1] == 'XTIME':
[365]	694	dtimeunits = var1.getncattr('description')
	695	tunits = dtimeunits.split(' ')[0]
	696	else:
	697	dtimeunits = var1.getncattr('units')
	698	tunits = dtimeunits.split(' ')[0]
	699
[1825]	700	dtime = (var1[1] - var1[0])*diag.timeunits_seconds(tunits)
[1908]	701	ncvar.insert_variable(ncobj, depvars[2], diagout/dtime, diagoutd, diagoutvd, \
	702	newnc)
[365]	703
[1909]	704	# fog_K84: Computation of fog and visibility following Kunkel, (1984) as QCLOUD, QICE
[1908]	705	elif diagn == 'fog_K84':
	706
	707	var0 = ncobj.variables[depvars[0]]
	708	var1 = ncobj.variables[depvars[1]]
	709
	710	dnamesvar = list(var0.dimensions)
	711	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	712
	713	diag1, diag2, diagoutd, diagoutvd = diag.Forcompute_fog_K84(var0, var1, \
	714	dnamesvar, dvnamesvar)
	715	# Removing the nonChecking variable-dimensions from the initial list
	716	varsadd = []
	717	diagoutvd = list(dvnames)
	718	for nonvd in NONchkvardims:
	719	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	720	varsadd.append(nonvd)
	721	ncvar.insert_variable(ncobj, 'fog', diag1, diagoutd, diagoutvd, newnc)
	722	ncvar.insert_variable(ncobj, 'fogvisblty', diag2, diagoutd, diagoutvd, newnc)
	723
[1909]	724	# fog_RUC: Computation of fog and visibility following Kunkel, (1984) as QVAPOR,
	725	# WRFt, WRFp or Q2, T2, PSFC
[1908]	726	elif diagn == 'fog_RUC':
	727
	728	var0 = ncobj.variables[depvars[0]]
[1909]	729	print gen.infmsg
	730	if depvars[1] == 'WRFt':
	731	print ' ' + main + ": computing '" + diagn + "' using 3D variables !!"
	732	var1 = WRFt
	733	var2 = WRFp
	734	else:
	735	print ' ' + main + ": computing '" + diagn + "' using 2D variables !!"
	736	var1 = ncobj.variables[depvars[1]]
	737	var2 = ncobj.variables[depvars[2]]
[1908]	738
	739	dnamesvar = list(var0.dimensions)
	740	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	741
[1909]	742	diag1, diag2, diagoutd, diagoutvd = diag.Forcompute_fog_RUC(var0, var1, var2,\
[1908]	743	dnamesvar, dvnamesvar)
	744	# Removing the nonChecking variable-dimensions from the initial list
	745	varsadd = []
	746	diagoutvd = list(dvnames)
	747	for nonvd in NONchkvardims:
	748	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	749	varsadd.append(nonvd)
	750	ncvar.insert_variable(ncobj, 'fog', diag1, diagoutd, diagoutvd, newnc)
	751	ncvar.insert_variable(ncobj, 'fogvisblty', diag2, diagoutd, diagoutvd, newnc)
	752
[1909]	753	# fog_FRAML50: Computation of fog and visibility following Gultepe, I. and
	754	# J.A. Milbrandt, 2010 as QVAPOR, WRFt, WRFp or Q2, T2, PSFC
	755	elif diagn == 'fog_FRAML50':
	756
	757	var0 = ncobj.variables[depvars[0]]
	758	print gen.infmsg
	759	if depvars[1] == 'WRFt':
	760	print ' ' + main + ": computing '" + diagn + "' using 3D variables !!"
	761	var1 = WRFt
	762	var2 = WRFp
	763	else:
	764	print ' ' + main + ": computing '" + diagn + "' using 2D variables !!"
	765	var1 = ncobj.variables[depvars[1]]
	766	var2 = ncobj.variables[depvars[2]]
	767
	768	dnamesvar = list(var0.dimensions)
	769	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	770
	771	diag1, diag2, diagoutd, diagoutvd = diag.Forcompute_fog_FRAML50(var0, var1, \
	772	var2, dnamesvar, dvnamesvar)
	773	# Removing the nonChecking variable-dimensions from the initial list
	774	varsadd = []
	775	diagoutvd = list(dvnames)
	776	for nonvd in NONchkvardims:
	777	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	778	varsadd.append(nonvd)
	779	ncvar.insert_variable(ncobj, 'fog', diag1, diagoutd, diagoutvd, newnc)
	780	ncvar.insert_variable(ncobj, 'fogvisblty', diag2, diagoutd, diagoutvd, newnc)
	781
[2643]	782	# front_R04: Computation of the presence of a front as defined by Rodrigues et al.
[2655]	783	# (2004) (tas, uas, vas, dxs, dys, 'params',[dtas],[dwss])
[2643]	784	# Rodrigues et al. 2004: Climatologia de frentes frias no litoral de Santa Catarina,
	785	# Rev. Bras. Geof. v.22 n.2 Sao Paulo maio/ago. DOI: 10.1590/S0102-261X2004000200004
	786	elif diagn == 'front_R04':
	787
	788	var0 = ncobj.variables[depvars[0]]
	789	var1 = ncobj.variables[depvars[1]]
	790	var2 = ncobj.variables[depvars[2]]
[2655]	791	if depvars[3] == 'WRFdx': var3 = WRFdx
	792	else: var3 = ncobj.variables[depvars[3]]
	793	if depvars[4] == 'WRFdy': var4 = WRFdy
	794	else: var4 = ncobj.variables[depvars[4]]
	795	par1 = np.float(depvars[5].split(':')[1])
	796	par2 = np.float(depvars[5].split(':')[2])
[2643]	797
	798	dnamesvar = list(var0.dimensions)
	799	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	800
[2655]	801	diag1, diag2, diag3, diag4, diagoutd, diagoutvd = diag.Forcompute_front_R04( \
	802	var0, var1, var2, var3, var4, par1, par2, dnamesvar, dvnamesvar)
[2643]	803
	804	# Removing the nonChecking variable-dimensions from the initial list
	805	varsadd = []
	806	diagoutvd = list(dvnames)
	807	for nonvd in NONchkvardims:
	808	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	809	varsadd.append(nonvd)
[2655]	810	ncvar.insert_variable(ncobj, 'front', diag1, diagoutd, diagoutvd, newnc, \
	811	fill=gen.fillValueF)
	812	ovar = newnc.variables['front']
	813	ovar.setncattr('description', 'front defintion after a generalization of ' + \
	814	'Rodrigues et al. 2004, Rev. Bras. Geof. v.22 n.2, ' + \
	815	'doi: 10.1590/S0102-261X2004000200004')
	816	ovar.setncattr('details1', 'N-S wind gradient replaced by gradient of uas, '+\
	817	'vas')
	818	ovar.setncattr('dtas', par1)
	819	ovar.setncattr('dwss', par2)
	820	newnc.sync()
	821	ncvar.insert_variable(ncobj, 'dt1', diag2, diagoutd, diagoutvd, newnc, \
	822	fill=gen.fillValueF)
	823	ovar = newnc.variables['dt1']
	824	vu = var0.getncattr('units')
	825	ncvar.set_attribute(ovar, 'units', vu)
	826	attrv = ovar.getncattr('standard_name')
	827	ncvar.set_attribute(ovar, 'standard_name', attrv + depvars[0])
	828	attrv = ovar.getncattr('long_name')
	829	ncvar.set_attribute(ovar, 'long_name', attrv + ' of ' + depvars[0])
	830	ncvar.insert_variable(ncobj, 'gradh', diag3, diagoutd, diagoutvd, newnc, \
	831	fill=gen.fillValueF)
	832	ovar = newnc.variables['gradh']
	833	vu = var1.getncattr('units')
	834	ncvar.set_attribute(ovar, 'units', vu+'m-1')
	835	attrv = ovar.getncattr('standard_name')
	836	ncvar.set_attribute(ovar, 'standard_name', attrv + 'wss')
	837	attrv = ovar.getncattr('long_name')
	838	ncvar.set_attribute(ovar, 'long_name', attrv + ' of wss')
	839	ncvar.insert_variable(ncobj, 'dt2', diag4, diagoutd, diagoutvd, newnc, \
	840	fill=gen.fillValueF)
	841	vu = var0.getncattr('units')
	842	ncvar.set_attribute(ovar, 'units', vu)
	843	attrv = ovar.getncattr('standard_name')
	844	ncvar.set_attribute(ovar, 'standard_name', attrv + depvars[0])
	845	attrv = ovar.getncattr('long_name')
	846	ncvar.set_attribute(ovar, 'long_name', attrv + ' of ' + depvars[0])
[2643]	847
[365]	848	# LMDZrh (pres, t, r)
[1758]	849	elif diagn == 'LMDZrh':
[365]	850
	851	var0 = ncobj.variables[depvars[0]][:]
	852	var1 = ncobj.variables[depvars[1]][:]
	853	var2 = ncobj.variables[depvars[2]][:]
	854
[1675]	855	diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnames,dvnames)
[1079]	856	ncvar.insert_variable(ncobj, 'hur', diagout, diagoutd, diagoutvd, newnc)
[365]	857
	858	# LMDZrhs (psol, t2m, q2m)
[1758]	859	elif diagn == 'LMDZrhs':
[365]	860
	861	var0 = ncobj.variables[depvars[0]][:]
	862	var1 = ncobj.variables[depvars[1]][:]
	863	var2 = ncobj.variables[depvars[2]][:]
	864
	865	dnamesvar = ncobj.variables[depvars[0]].dimensions
	866	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	867
[1675]	868	diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
[365]	869
[1079]	870	ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
[365]	871
[2215]	872	# range_faces: LON, LAT, HGT, WRFdxdy, 'face:['WE'/'SN']:[dsfilt]:[dsnewrange]:[hvalleyrange]'
[2208]	873	elif diagn == 'range_faces':
	874
	875	var0 = ncobj.variables[depvars[0]][:]
	876	var1 = ncobj.variables[depvars[1]][:]
	877	var2 = ncobj.variables[depvars[2]][:]
[2215]	878	face = depvars[4].split(':')[1]
	879	dsfilt = np.float(depvars[4].split(':')[2])
	880	dsnewrange = np.float(depvars[4].split(':')[3])
	881	hvalleyrange = np.float(depvars[4].split(':')[4])
[2208]	882
[2222]	883	dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
[2212]	884	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
[2209]	885	lon, lat = gen.lonlat2D(var0, var1)
	886	if len(var2.shape) == 3:
	887	print warnmsg
	888	print ' ' + diagn + ": shapping to 2D variable '" + depvars[2] + "' !!"
	889	hgt = var2[0,:,:]
[2212]	890	dnamesvar.pop(0)
	891	dvnamesvar.pop(0)
[2209]	892	else:
	893	hgt = var2[:]
	894
[2213]	895	orogmax, ptorogmax, dhgt, peaks, valleys, origfaces, diagout, diagoutd, \
[2223]	896	diagoutvd, rng, rngorogmax, ptrngorogmax= diag.Forcompute_range_faces(lon, \
[2260]	897	lat, hgt, WRFdx, WRFdy, WRFds, face, dsfilt, dsnewrange, hvalleyrange, \
	898	dnamesvar, dvnamesvar)
[2208]	899
	900	# Removing the nonChecking variable-dimensions from the initial list
	901	varsadd = []
	902	diagoutvd = list(dvnames)
	903	for nonvd in NONchkvardims:
	904	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	905	varsadd.append(nonvd)
[2212]	906
[2215]	907	ncvar.insert_variable(ncobj, 'dx', WRFdx, diagoutd, diagoutvd, newnc)
	908	ncvar.insert_variable(ncobj, 'dy', WRFdy, diagoutd, diagoutvd, newnc)
	909	ncvar.insert_variable(ncobj, 'ds', WRFds, diagoutd, diagoutvd, newnc)
	910
[2213]	911	# adding variables to output file
	912	if face == 'WE': axis = 'lon'
	913	elif face == 'SN': axis = 'lat'
	914
[2223]	915	ncvar.insert_variable(ncobj, 'range', rng, diagoutd, diagoutvd, newnc, \
	916	fill=gen.fillValueI)
	917	ovar = newnc.variables['range']
	918	ncvar.set_attribute(ovar, 'deriv', axis)
	919
[2214]	920	ncvar.insert_variable(ncobj, 'orogmax', rngorogmax, diagoutd, diagoutvd, \
[2222]	921	newnc, fill=gen.fillValueF)
[2214]	922	newnc.renameVariable('orogmax', 'rangeorogmax')
	923	ovar = newnc.variables['rangeorogmax']
	924	ncvar.set_attribute(ovar, 'deriv', axis)
	925	stdn = ovar.standard_name
	926	ncvar.set_attribute(ovar, 'standard_name', 'range_' + stdn)
	927	Ln = ovar.long_name
	928	ncvar.set_attribute(ovar, 'long_name', 'range ' + stdn)
	929
	930	ncvar.insert_variable(ncobj, 'ptorogmax', ptrngorogmax, diagoutd, diagoutvd, \
	931	newnc)
	932	newnc.renameVariable('ptorogmax', 'rangeptorogmax')
	933	ovar = newnc.variables['rangeptorogmax']
	934	ncvar.set_attribute(ovar, 'deriv', axis)
	935	stdn = ovar.standard_name
	936	ncvar.set_attribute(ovar, 'standard_name', 'range_' + stdn)
	937	Ln = ovar.long_name
	938	ncvar.set_attribute(ovar, 'long_name', 'range ' + stdn)
	939
[2213]	940	ncvar.insert_variable(ncobj, 'orogmax', orogmax, diagoutd, diagoutvd, \
	941	newnc)
	942	ovar = newnc.variables['orogmax']
	943	ncvar.set_attribute(ovar, 'deriv', axis)
	944
	945	ncvar.insert_variable(ncobj, 'ptorogmax', ptorogmax, diagoutd, diagoutvd, \
	946	newnc)
	947	ovar = newnc.variables['ptorogmax']
	948	ncvar.set_attribute(ovar, 'deriv', axis)
	949
[2215]	950	ncvar.insert_variable(ncobj, 'orogderiv', dhgt, diagoutd, diagoutvd, newnc)
	951	ovar = newnc.variables['orogderiv']
[2212]	952	ncvar.set_attribute(ovar, 'deriv', axis)
[2208]	953
[2212]	954	ncvar.insert_variable(ncobj, 'peak', peaks, diagoutd, diagoutvd, newnc)
	955	ncvar.insert_variable(ncobj, 'valley', valleys, diagoutd, diagoutvd, newnc)
	956
	957	ncvar.insert_variable(ncobj, 'rangefaces', diagout, diagoutd, diagoutvd, \
	958	newnc)
	959	newnc.renameVariable('rangefaces', 'rangefacesfilt')
	960	ovar = newnc.variables['rangefacesfilt']
	961	ncvar.set_attribute(ovar, 'face', face)
[2215]	962	ncvar.set_attributek(ovar, 'dist_filter', dsfilt, 'F')
[2212]	963
	964	ncvar.insert_variable(ncobj, 'rangefaces', origfaces, diagoutd, diagoutvd, \
[2215]	965	newnc, fill=gen.fillValueI)
[2212]	966	ovar = newnc.variables['rangefaces']
	967	ncvar.set_attribute(ovar, 'face', face)
[2215]	968	ncvar.set_attributek(ovar, 'dist_newrange', dsnewrange, 'F')
	969	ncvar.set_attributek(ovar, 'h_valley_newrange', hvalleyrange, 'F')
[2212]	970
[2277]	971	# cell_bnds: grid cell bounds from lon, lat from a reglar lon/lat projection as
	972	# intersection of their related parallels and meridians
	973	elif diagn == 'reglonlatbnds':
	974
	975	var00 = ncobj.variables[depvars[0]][:]
	976	var01 = ncobj.variables[depvars[1]][:]
	977
	978	var0, var1 = gen.lonlat2D(var00,var01)
	979
	980	dnamesvar = []
	981	dnamesvar.append('bnds')
	982	if (len(var00.shape) == 3):
	983	dnamesvar.append(ncobj.variables[depvars[0]].dimensions[1])
	984	dnamesvar.append(ncobj.variables[depvars[0]].dimensions[2])
	985	elif (len(var00.shape) == 2):
	986	dnamesvar.append(ncobj.variables[depvars[0]].dimensions[0])
	987	dnamesvar.append(ncobj.variables[depvars[0]].dimensions[1])
	988	elif (len(var00.shape) == 1):
	989	dnamesvar.append(ncobj.variables[depvars[0]].dimensions[0])
	990	dnamesvar.append(ncobj.variables[depvars[1]].dimensions[0])
	991	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	992
	993	cellbndsx, cellbndsy, diagoutd, diagoutvd = diag.Forcompute_cellbndsreg(var0,\
	994	var1, dnamesvar, dvnamesvar)
	995
	996	# Removing the nonChecking variable-dimensions from the initial list
	997	varsadd = []
	998	diagoutvd = list(dvnames)
	999	for nonvd in NONchkvardims:
	1000	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1001	varsadd.append(nonvd)
	1002	# creation of bounds dimension
	1003	newdim = newnc.createDimension('bnds', 4)
	1004
	1005	ncvar.insert_variable(ncobj, 'lon_bnds', cellbndsx, diagoutd, diagoutvd, newnc)
	1006	ncvar.insert_variable(ncobj, 'lat_bnds', cellbndsy, diagoutd, diagoutvd, newnc)
	1007
[2390]	1008	# tws: Bet Wulb temperature following Stull 2011 (tas, hurs)
	1009	elif diagn == 'tws':
[2387]	1010
	1011	var0 = ncobj.variables[depvars[0]][:]
	1012	var1 = ncobj.variables[depvars[1]][:]
	1013
[2390]	1014	dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
	1015	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1016	diagoutd = dnames
	1017	diagoutvd = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1018
	1019	diagout = diag.var_tws_S11(var0,var1)
[2387]	1020	ncvar.insert_variable(ncobj, 'tws', diagout, diagoutd, diagoutvd, newnc)
	1021
[2274]	1022	# cell_bnds: grid cell bounds from XLONG_U, XLAT_U, XLONG_V, XLAT_V as intersection
	1023	# of their related parallels and meridians
	1024	elif diagn == 'WRFbnds':
	1025
	1026	var0 = ncobj.variables[depvars[0]][0,:,:]
	1027	var1 = ncobj.variables[depvars[1]][0,:,:]
	1028	var2 = ncobj.variables[depvars[2]][0,:,:]
	1029	var3 = ncobj.variables[depvars[3]][0,:,:]
	1030
	1031	dnamesvar = []
[2276]	1032	dnamesvar.append('bnds')
	1033	dnamesvar.append(ncobj.variables[depvars[0]].dimensions[1])
[2274]	1034	dnamesvar.append(ncobj.variables[depvars[2]].dimensions[2])
	1035	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1036
	1037	cellbndsx, cellbndsy, diagoutd, diagoutvd = diag.Forcompute_cellbnds(var0, \
	1038	var1, var2, var3, dnamesvar, dvnamesvar)
	1039
	1040	# Removing the nonChecking variable-dimensions from the initial list
	1041	varsadd = []
	1042	diagoutvd = list(dvnames)
	1043	for nonvd in NONchkvardims:
	1044	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1045	varsadd.append(nonvd)
	1046	# creation of bounds dimension
	1047	newdim = newnc.createDimension('bnds', 4)
	1048
	1049	ncvar.insert_variable(ncobj, 'lon_bnds', cellbndsx, diagoutd, diagoutvd, newnc)
[2285]	1050	newnc.sync()
[2274]	1051	ncvar.insert_variable(ncobj, 'lat_bnds', cellbndsy, diagoutd, diagoutvd, newnc)
[2285]	1052	newnc.sync()
[2274]	1053
[2346]	1054	if newnc.variables.has_key('XLONG'):
	1055	ovar = newnc.variables['XLONG']
	1056	ovar.setncattr('bounds', 'lon_bnds lat_bnds')
	1057	ovar = newnc.variables['XLAT']
	1058	ovar.setncattr('bounds', 'lon_bnds lat_bnds')
	1059	elif newnc.variables.has_key('XLONG_M'):
	1060	ovar = newnc.variables['XLONG_M']
	1061	ovar.setncattr('bounds', 'lon_bnds lat_bnds')
	1062	ovar = newnc.variables['XLAT_M']
	1063	ovar.setncattr('bounds', 'lon_bnds lat_bnds')
	1064	else:
	1065	print errormsg
	1066	print ' ' + fname + ": error computing diagnostic '" + diagn + "' !!"
	1067	print " neither: 'XLONG', 'XLONG_M' have been found"
	1068	quit(-1)
[2283]	1069
[1762]	1070	# mrso: total soil moisture SMOIS, DZS
	1071	elif diagn == 'WRFmrso':
	1072
	1073	var0 = ncobj.variables[depvars[0]][:]
	1074	var10 = ncobj.variables[depvars[1]][:]
	1075	dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
	1076	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1077
	1078	var1 = var0.copy()*0.
	1079	var2 = var0.copy()*0.+1.
	1080	# Must be a better way....
	1081	for j in range(var0.shape[2]):
	1082	for i in range(var0.shape[3]):
	1083	var1[:,:,j,i] = var10
	1084
	1085	diagout, diagoutd, diagoutvd = diag.Forcompute_zint(var0, var1, var2, \
	1086	dnamesvar, dvnamesvar)
	1087
	1088	# Removing the nonChecking variable-dimensions from the initial list
	1089	varsadd = []
	1090	diagoutvd = list(dvnames)
	1091	for nonvd in NONchkvardims:
	1092	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1093	varsadd.append(nonvd)
	1094	ncvar.insert_variable(ncobj, 'mrso', diagout, diagoutd, diagoutvd, newnc)
	1095
[2206]	1096	# mrsos: First layer soil moisture SMOIS, DZS
	1097	elif diagn == 'WRFmrsos':
	1098
	1099	var0 = ncobj.variables[depvars[0]][:]
	1100	var1 = ncobj.variables[depvars[1]][:]
	1101	diagoutd = list(ncobj.variables[depvars[0]].dimensions)
	1102	diagoutvd = ncvar.var_dim_dimv(diagoutd,dnames,dvnames)
	1103
	1104	diagoutd.pop(1)
	1105	diagoutvd.pop(1)
	1106
	1107	diagout= np.zeros((var0.shape[0],var0.shape[2],var0.shape[3]), dtype=np.float)
	1108
	1109	# Must be a better way....
	1110	for j in range(var0.shape[2]):
	1111	for i in range(var0.shape[3]):
	1112	diagout[:,j,i] = var0[:,0,j,i]*var1[:,0]
	1113
	1114	# Removing the nonChecking variable-dimensions from the initial list
	1115	varsadd = []
	1116	diagoutvd = list(dvnames)
	1117	for nonvd in NONchkvardims:
	1118	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1119	varsadd.append(nonvd)
	1120	ncvar.insert_variable(ncobj, 'mrsos', diagout, diagoutd, diagoutvd, newnc)
	1121
[365]	1122	# mslp: mean sea level pressure (pres, psfc, terrain, temp, qv)
[1758]	1123	elif diagn == 'mslp' or diagn == 'WRFmslp':
[365]	1124
	1125	var1 = ncobj.variables[depvars[1]][:]
	1126	var2 = ncobj.variables[depvars[2]][:]
	1127	var4 = ncobj.variables[depvars[4]][:]
	1128
[1758]	1129	if diagn == 'WRFmslp':
[365]	1130	var0 = WRFp
	1131	var3 = WRFt
	1132	dnamesvar = ncobj.variables['P'].dimensions
	1133	else:
	1134	var0 = ncobj.variables[depvars[0]][:]
	1135	var3 = ncobj.variables[depvars[3]][:]
	1136	dnamesvar = ncobj.variables[depvars[0]].dimensions
	1137
	1138	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1139
[1675]	1140	diagout, diagoutd, diagoutvd = diag.compute_mslp(var0, var1, var2, var3, var4, \
[365]	1141	dnamesvar, dvnamesvar)
	1142
[1581]	1143	# Removing the nonChecking variable-dimensions from the initial list
	1144	varsadd = []
	1145	diagoutvd = list(dvnames)
	1146	for nonvd in NONchkvardims:
	1147	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1148	varsadd.append(nonvd)
[365]	1149	ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
	1150
[2387]	1151	# WRFtws: Bet Wulb temperature following Stull 2011 (PSFC, T2, Q2)
	1152	elif diagn == 'WRFtws':
	1153
	1154	var0 = ncobj.variables[depvars[0]][:]
	1155	var1 = ncobj.variables[depvars[1]][:]
	1156	var2 = ncobj.variables[depvars[2]][:]
	1157
	1158	dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
	1159	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1160
	1161	hurs, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
	1162
	1163	diagout = diag.var_tws_S11(var1, hurs)
	1164	# Removing the nonChecking variable-dimensions from the initial list
	1165	varsadd = []
	1166	diagoutvd = list(dvnames)
	1167	for nonvd in NONchkvardims:
	1168	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1169	varsadd.append(nonvd)
	1170
	1171	ncvar.insert_variable(ncobj, 'tws', diagout, diagoutd, diagoutvd, newnc)
	1172
[642]	1173	# OMEGAw (omega, p, t) from NCL formulation (https://www.ncl.ucar.edu/Document/Functions/Contributed/omega_to_w.shtml)
[1758]	1174	elif diagn == 'OMEGAw':
[642]	1175
	1176	var0 = ncobj.variables[depvars[0]][:]
	1177	var1 = ncobj.variables[depvars[1]][:]
[643]	1178	var2 = ncobj.variables[depvars[2]][:]
[642]	1179
	1180	dnamesvar = ncobj.variables[depvars[0]].dimensions
	1181	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1182
[1675]	1183	diagout, diagoutd, diagoutvd = diag.compute_OMEGAw(var0,var1,var2,dnamesvar,dvnamesvar)
[642]	1184
	1185	ncvar.insert_variable(ncobj, 'wa', diagout, diagoutd, diagoutvd, newnc)
	1186
[2095]	1187	# raintot: instantaneous total precipitation from WRF as (RAINC + RAINC + RAINSH) / dTime
[1758]	1188	elif diagn == 'RAINTOT':
[365]	1189
	1190	var0 = ncobj.variables[depvars[0]]
	1191	var1 = ncobj.variables[depvars[1]]
[2095]	1192	var2 = ncobj.variables[depvars[2]]
	1193
	1194	if depvars[3] != 'WRFtime':
	1195	var3 = ncobj.variables[depvars[3]]
[654]	1196	else:
[2095]	1197	var3 = np.arange(var0.shape[0], dtype=int)
[365]	1198
[2095]	1199	var = var0[:] + var1[:] + var2[:]
[365]	1200
	1201	dnamesvar = var0.dimensions
	1202	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1203
[1675]	1204	diagout, diagoutd, diagoutvd = diag.compute_deaccum(var,dnamesvar,dvnamesvar)
[365]	1205
	1206	# Transforming to a flux
[2095]	1207	if var3.shape[0] > 1:
	1208	if depvars[3] != 'WRFtime':
	1209	dtimeunits = var3.getncattr('units')
[600]	1210	tunits = dtimeunits.split(' ')[0]
	1211
[2095]	1212	dtime = (var3[1] - var3[0])*diag.timeunits_seconds(tunits)
[600]	1213	else:
[2095]	1214	var3 = ncobj.variables['Times']
	1215	time1 = var3[0,:]
	1216	time2 = var3[1,:]
[600]	1217	tmf1 = ''
	1218	tmf2 = ''
	1219	for ic in range(len(time1)):
	1220	tmf1 = tmf1 + time1[ic]
	1221	tmf2 = tmf2 + time2[ic]
[654]	1222	dtdate1 = dtime.datetime.strptime(tmf1,"%Y-%m-%d_%H:%M:%S")
	1223	dtdate2 = dtime.datetime.strptime(tmf2,"%Y-%m-%d_%H:%M:%S")
[600]	1224	diffdate12 = dtdate2 - dtdate1
	1225	dtime = diffdate12.total_seconds()
	1226	print 'dtime:',dtime
[442]	1227	else:
[600]	1228	print warnmsg
[1645]	1229	print ' ' + main + ": only 1 time-step for '" + diag + "' !!"
[600]	1230	print ' leaving a zero value!'
[1646]	1231	diagout = var0[:]*0.
[600]	1232	dtime=1.
[442]	1233
[1644]	1234	# Removing the nonChecking variable-dimensions from the initial list
	1235	varsadd = []
	1236	for nonvd in NONchkvardims:
	1237	if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
	1238	varsadd.append(nonvd)
	1239
[365]	1240	ncvar.insert_variable(ncobj, 'pr', diagout/dtime, diagoutd, diagoutvd, newnc)
	1241
[2140]	1242	# timemax ([varname], time). When a given variable [varname] got its maximum
	1243	elif diagn == 'timemax':
	1244
	1245	var0 = ncobj.variables[depvars[0]][:]
	1246	var1 = ncobj.variables[depvars[1]][:]
	1247
	1248	otime = ncobj.variables[depvars[1]]
	1249
	1250	dnamesvar = ncobj.variables[depvars[0]].dimensions
	1251	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1252
	1253	diagout, diagoutd, diagoutvd = diag.var_timemax(var0, var1, dnames, \
	1254	dvnames)
	1255
	1256	ncvar.insert_variable(ncobj, 'timemax', diagout, diagoutd, diagoutvd, newnc, \
	1257	fill=gen.fillValueF)
	1258	# Getting the right units
	1259	ovar = newnc.variables['timemax']
	1260	if gen.searchInlist(otime.ncattrs(), 'units'):
	1261	tunits = otime.getncattr('units')
	1262	ncvar.set_attribute(ovar, 'units', tunits)
	1263	newnc.sync()
	1264	ncvar.set_attribute(ovar, 'variable', depvars[0])
	1265
[2138]	1266	# timeoverthres ([varname], time, [value], [CFvarn]). When a given variable [varname]
	1267	# overpass a given [value]. Being [CFvarn] the name of the diagnostics in
	1268	# variables_values.dat
	1269	elif diagn == 'timeoverthres':
	1270
	1271	var0 = ncobj.variables[depvars[0]][:]
	1272	var1 = ncobj.variables[depvars[1]][:]
	1273	var2 = np.float(depvars[2])
	1274	var3 = depvars[3]
	1275
	1276	otime = ncobj.variables[depvars[1]]
	1277
	1278	dnamesvar = ncobj.variables[depvars[0]].dimensions
	1279	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1280
	1281	diagout, diagoutd, diagoutvd = diag.var_timeoverthres(var0, var1, var2, \
	1282	dnames, dvnames)
	1283
	1284	ncvar.insert_variable(ncobj, var3, diagout, diagoutd, diagoutvd, newnc, \
	1285	fill=gen.fillValueF)
	1286	# Getting the right units
	1287	ovar = newnc.variables[var3]
	1288	if gen.searchInlist(otime.ncattrs(), 'units'):
	1289	tunits = otime.getncattr('units')
	1290	ncvar.set_attribute(ovar, 'units', tunits)
	1291	newnc.sync()
	1292	ncvar.set_attribute(ovar, 'overpassed_threshold', var2)
	1293
[612]	1294	# rhs (psfc, t, q) from TimeSeries files
[1758]	1295	elif diagn == 'TSrhs':
[612]	1296
	1297	p0=100000.
	1298	var0 = ncobj.variables[depvars[0]][:]
	1299	var1 = (ncobj.variables[depvars[1]][:])(var0/p0)*(2./7.)
	1300	var2 = ncobj.variables[depvars[2]][:]
	1301
	1302	dnamesvar = ncobj.variables[depvars[0]].dimensions
	1303	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1304
[1675]	1305	diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
[612]	1306
[1079]	1307	ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
[612]	1308
[2390]	1309	# rhs (psfc, t, q) from tas, tds
	1310	elif diagn == 'rhs_tas_tds':
	1311
	1312	var0 = ncobj.variables[depvars[0]][:]
	1313	var1 = ncobj.variables[depvars[1]][:]
	1314
	1315	dnamesvar = ncobj.variables[depvars[0]].dimensions
	1316	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1317
	1318	diagout, diagoutd, diagoutvd = diag.var_hur_tas_tds(var0,var1,dnamesvar, \
	1319	dvnamesvar)
	1320
	1321	ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
	1322
[1762]	1323	# slw: total soil liquid water SH2O, DZS
	1324	elif diagn == 'WRFslw':
	1325
	1326	var0 = ncobj.variables[depvars[0]][:]
	1327	var10 = ncobj.variables[depvars[1]][:]
	1328	dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
	1329	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1330
	1331	var1 = var0.copy()*0.
	1332	var2 = var0.copy()*0.+1.
	1333	# Must be a better way....
	1334	for j in range(var0.shape[2]):
	1335	for i in range(var0.shape[3]):
	1336	var1[:,:,j,i] = var10
	1337
	1338	diagout, diagoutd, diagoutvd = diag.Forcompute_zint(var0, var1, var2, \
	1339	dnamesvar, dvnamesvar)
	1340
	1341	# Removing the nonChecking variable-dimensions from the initial list
	1342	varsadd = []
	1343	diagoutvd = list(dvnames)
	1344	for nonvd in NONchkvardims:
	1345	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1346	varsadd.append(nonvd)
	1347	ncvar.insert_variable(ncobj, 'slw', diagout, diagoutd, diagoutvd, newnc)
	1348
[612]	1349	# td (psfc, t, q) from TimeSeries files
[1758]	1350	elif diagn == 'TStd' or diagn == 'td':
[612]	1351
	1352	var0 = ncobj.variables[depvars[0]][:]
	1353	var1 = ncobj.variables[depvars[1]][:] - 273.15
	1354	var2 = ncobj.variables[depvars[2]][:]
	1355
	1356	dnamesvar = ncobj.variables[depvars[0]].dimensions
	1357	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1358
[1675]	1359	diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
[612]	1360
[1966]	1361	ncvar.insert_variable(ncobj, 'tdas', diagout, diagoutd, diagoutvd, newnc)
[612]	1362
	1363	# td (psfc, t, q) from TimeSeries files
[1758]	1364	elif diagn == 'TStdC' or diagn == 'tdC':
[612]	1365
	1366	var0 = ncobj.variables[depvars[0]][:]
	1367	# Temperature is already in degrees Celsius
	1368	var1 = ncobj.variables[depvars[1]][:]
	1369	var2 = ncobj.variables[depvars[2]][:]
	1370
	1371	dnamesvar = ncobj.variables[depvars[0]].dimensions
	1372	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1373
[1675]	1374	diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
[612]	1375
[1999]	1376	# Removing the nonChecking variable-dimensions from the initial list
	1377	varsadd = []
	1378	for nonvd in NONchkvardims:
	1379	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1380	varsadd.append(nonvd)
	1381
[1966]	1382	ncvar.insert_variable(ncobj, 'tdas', diagout, diagoutd, diagoutvd, newnc)
[612]	1383
	1384	# wds (u, v)
[1758]	1385	elif diagn == 'TSwds' or diagn == 'wds' :
[612]	1386
	1387	var0 = ncobj.variables[depvars[0]][:]
	1388	var1 = ncobj.variables[depvars[1]][:]
	1389
	1390	dnamesvar = ncobj.variables[depvars[0]].dimensions
	1391	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1392
[1675]	1393	diagout, diagoutd, diagoutvd = diag.compute_wds(var0,var1,dnamesvar,dvnamesvar)
[612]	1394
[1999]	1395	# Removing the nonChecking variable-dimensions from the initial list
	1396	varsadd = []
	1397	for nonvd in NONchkvardims:
	1398	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1399	varsadd.append(nonvd)
	1400
[612]	1401	ncvar.insert_variable(ncobj, 'wds', diagout, diagoutd, diagoutvd, newnc)
	1402
	1403	# wss (u, v)
[1758]	1404	elif diagn == 'TSwss' or diagn == 'wss':
[612]	1405
	1406	var0 = ncobj.variables[depvars[0]][:]
	1407	var1 = ncobj.variables[depvars[1]][:]
	1408
	1409	dnamesvar = ncobj.variables[depvars[0]].dimensions
	1410	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1411
[1675]	1412	diagout, diagoutd, diagoutvd = diag.compute_wss(var0,var1,dnamesvar,dvnamesvar)
[612]	1413
[1999]	1414	# Removing the nonChecking variable-dimensions from the initial list
	1415	varsadd = []
	1416	for nonvd in NONchkvardims:
	1417	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1418	varsadd.append(nonvd)
	1419
[612]	1420	ncvar.insert_variable(ncobj, 'wss', diagout, diagoutd, diagoutvd, newnc)
	1421
[365]	1422	# turbulence (var)
[1758]	1423	elif diagn == 'turbulence':
[365]	1424
	1425	var0 = ncobj.variables[depvars][:]
	1426
	1427	dnamesvar = list(ncobj.variables[depvars].dimensions)
	1428	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1429
[1675]	1430	diagout, diagoutd, diagoutvd = diag.compute_turbulence(var0,dnamesvar,dvnamesvar)
[959]	1431	valsvar = gen.variables_values(depvars)
[365]	1432
[959]	1433	newvarn = depvars + 'turb'
	1434	ncvar.insert_variable(ncobj, newvarn, diagout, diagoutd,
[365]	1435	diagoutvd, newnc)
[959]	1436	varobj = newnc.variables[newvarn]
[365]	1437	attrv = varobj.long_name
	1438	attr = varobj.delncattr('long_name')
	1439	newattr = ncvar.set_attribute(varobj, 'long_name', attrv + \
	1440	" Taylor decomposition turbulence term")
	1441
[1927]	1442	# ua va from ws wd (deg)
	1443	elif diagn == 'uavaFROMwswd':
	1444
	1445	var0 = ncobj.variables[depvars[0]][:]
	1446	var1 = ncobj.variables[depvars[1]][:]
	1447
	1448	ua = var0np.cos(var1np.pi/180.)
	1449	va = var0np.sin(var1np.pi/180.)
	1450
	1451	dnamesvar = ncobj.variables[depvars[0]].dimensions
	1452	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1453
	1454	ncvar.insert_variable(ncobj, 'ua', ua, dnamesvar, dvnamesvar, newnc)
	1455	ncvar.insert_variable(ncobj, 'va', va, dnamesvar, dvnamesvar, newnc)
	1456
[2033]	1457	# ua va from obs ws wd (deg)
	1458	elif diagn == 'uavaFROMobswswd':
	1459
	1460	var0 = ncobj.variables[depvars[0]][:]
	1461	var1 = ncobj.variables[depvars[1]][:]
[1927]	1462
[2033]	1463	ua = var0np.cos((var1+180.)np.pi/180.)
	1464	va = var0np.sin((var1+180.)np.pi/180.)
	1465
	1466	dnamesvar = ncobj.variables[depvars[0]].dimensions
	1467	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1468
	1469	ncvar.insert_variable(ncobj, 'ua', ua, dnamesvar, dvnamesvar, newnc)
	1470	ncvar.insert_variable(ncobj, 'va', va, dnamesvar, dvnamesvar, newnc)
	1471
[390]	1472	# WRFbils fom WRF as HFX + LH
[1758]	1473	elif diagn == 'WRFbils':
[390]	1474
	1475	var0 = ncobj.variables[depvars[0]][:]
	1476	var1 = ncobj.variables[depvars[1]][:]
	1477
	1478	diagout = var0 + var1
[867]	1479	dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
	1480	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
[390]	1481
[867]	1482	ncvar.insert_variable(ncobj, 'bils', diagout, dnamesvar, dvnamesvar, newnc)
[390]	1483
[1759]	1484	# WRFcape_afwa CAPE, CIN, ZLFC, PLFC, LI following WRF 'phys/module_diaf_afwa.F'
	1485	# methodology as WRFt, WRFrh, WRFp, WRFgeop, HGT
	1486	elif diagn == 'WRFcape_afwa':
	1487	var0 = WRFt
	1488	var1 = WRFrh
	1489	var2 = WRFp
	1490	dz = WRFgeop.shape[1]
	1491	# de-staggering
[1833]	1492	var3 = 0.5*(WRFgeop[:,0:dz-1,:,:]+WRFgeop[:,1:dz,:,:])/9.8
[1759]	1493	var4 = ncobj.variables[depvars[4]][0,:,:]
	1494
	1495	dnamesvar = list(ncobj.variables['T'].dimensions)
	1496	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1497
	1498	diagout = np.zeros(var0.shape, dtype=np.float)
	1499	diagout1, diagout2, diagout3, diagout4, diagout5, diagoutd, diagoutvd = \
	1500	diag.Forcompute_cape_afwa(var0, var1, var2, var3, var4, 3, dnamesvar, \
	1501	dvnamesvar)
	1502
	1503	# Removing the nonChecking variable-dimensions from the initial list
	1504	varsadd = []
	1505	for nonvd in NONchkvardims:
	1506	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1507	varsadd.append(nonvd)
	1508
	1509	ncvar.insert_variable(ncobj, 'cape', diagout1, diagoutd, diagoutvd, newnc)
	1510	ncvar.insert_variable(ncobj, 'cin', diagout2, diagoutd, diagoutvd, newnc)
	1511	ncvar.insert_variable(ncobj, 'zlfc', diagout3, diagoutd, diagoutvd, newnc)
	1512	ncvar.insert_variable(ncobj, 'plfc', diagout4, diagoutd, diagoutvd, newnc)
	1513	ncvar.insert_variable(ncobj, 'li', diagout5, diagoutd, diagoutvd, newnc)
	1514
[1581]	1515	# WRFclivi WRF water vapour path WRFdens, QICE, QGRAUPEL, QHAIL
[1758]	1516	elif diagn == 'WRFclivi':
[1581]	1517
	1518	var0 = WRFdens
	1519	qtot = ncobj.variables[depvars[1]]
	1520	qtotv = qtot[:]
	1521	Nspecies = len(depvars) - 2
	1522	for iv in range(Nspecies):
[1585]	1523	if ncobj.variables.has_key(depvars[iv+2]):
	1524	var1 = ncobj.variables[depvars[iv+2]][:]
	1525	qtotv = qtotv + var1
[1581]	1526
	1527	dnamesvar = list(qtot.dimensions)
	1528	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1529
[1675]	1530	diagout, diagoutd, diagoutvd = diag.compute_clivi(var0, qtotv, dnamesvar,dvnamesvar)
[1581]	1531
	1532	# Removing the nonChecking variable-dimensions from the initial list
	1533	varsadd = []
	1534	diagoutvd = list(dvnames)
	1535	for nonvd in NONchkvardims:
	1536	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1537	varsadd.append(nonvd)
	1538	ncvar.insert_variable(ncobj, 'clivi', diagout, diagoutd, diagoutvd, newnc)
	1539
[1803]	1540	# WRFclwvi WRF water cloud-condensed path WRFdens, QCLOUD, QICE, QGRAUPEL, QHAIL
	1541	elif diagn == 'WRFclwvi':
[1581]	1542
	1543	var0 = WRFdens
	1544	qtot = ncobj.variables[depvars[1]]
	1545	qtotv = ncobj.variables[depvars[1]]
	1546	Nspecies = len(depvars) - 2
	1547	for iv in range(Nspecies):
[1585]	1548	if ncobj.variables.has_key(depvars[iv+2]):
	1549	var1 = ncobj.variables[depvars[iv+2]]
	1550	qtotv = qtotv + var1[:]
[1581]	1551
	1552	dnamesvar = list(qtot.dimensions)
	1553	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1554
[1675]	1555	diagout, diagoutd, diagoutvd = diag.compute_clwvl(var0, qtotv, dnamesvar,dvnamesvar)
[1581]	1556
	1557	# Removing the nonChecking variable-dimensions from the initial list
	1558	varsadd = []
	1559	diagoutvd = list(dvnames)
	1560	for nonvd in NONchkvardims:
	1561	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1562	varsadd.append(nonvd)
[1803]	1563	ncvar.insert_variable(ncobj, 'clwvi', diagout, diagoutd, diagoutvd, newnc)
[1581]	1564
[1809]	1565	# WRF_denszint WRF vertical integration as WRFdens, sum(Q[water species1], ..., Q[water speciesN]), varn=[varN]
	1566	elif diagn == 'WRF_denszint':
	1567
	1568	var0 = WRFdens
	1569	varn = depvars[1].split('=')[1]
	1570	qtot = ncobj.variables[depvars[2]]
	1571	qtotv = ncobj.variables[depvars[2]]
	1572	Nspecies = len(depvars) - 2
	1573	for iv in range(Nspecies):
	1574	if ncobj.variables.has_key(depvars[iv+2]):
	1575	var1 = ncobj.variables[depvars[iv+2]]
	1576	qtotv = qtotv + var1[:]
	1577
	1578	dnamesvar = list(qtot.dimensions)
	1579	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1580
	1581	diagout, diagoutd, diagoutvd = diag.compute_clwvl(var0, qtotv, dnamesvar,dvnamesvar)
	1582
	1583	# Removing the nonChecking variable-dimensions from the initial list
	1584	varsadd = []
	1585	diagoutvd = list(dvnames)
	1586	for nonvd in NONchkvardims:
	1587	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1588	varsadd.append(nonvd)
	1589	ncvar.insert_variable(ncobj, varn, diagout, diagoutd, diagoutvd, newnc)
	1590
[654]	1591	# WRFgeop geopotential from WRF as PH + PHB
[1758]	1592	elif diagn == 'WRFgeop':
[1382]	1593	var0 = ncobj.variables[depvars[0]][:]
	1594	var1 = ncobj.variables[depvars[1]][:]
[654]	1595
[1382]	1596	# de-staggering geopotential
	1597	diagout0 = var0 + var1
	1598	dt = diagout0.shape[0]
	1599	dz = diagout0.shape[1]
	1600	dy = diagout0.shape[2]
	1601	dx = diagout0.shape[3]
	1602
	1603	diagout = np.zeros((dt,dz-1,dy,dx), dtype=np.float)
	1604	diagout = 0.5*(diagout0[:,1:dz,:,:]+diagout0[:,0:dz-1,:,:])
	1605
	1606	# Removing the nonChecking variable-dimensions from the initial list
	1607	varsadd = []
[1389]	1608	diagoutvd = list(dvnames)
[1382]	1609	for nonvd in NONchkvardims:
[1389]	1610	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
[1382]	1611	varsadd.append(nonvd)
	1612
[1389]	1613	ncvar.insert_variable(ncobj, 'zg', diagout, dnames, diagoutvd, newnc)
[654]	1614
[1804]	1615	# WRFpotevap_orPM potential evapotranspiration following Penman-Monteith formulation
[1833]	1616	# implemented in ORCHIDEE (in src_sechiba/enerbil.f90) as: WRFdens, UST, U10, V10, T2, PSFC, QVAPOR
[1804]	1617	elif diagn == 'WRFpotevap_orPM':
	1618	var0 = WRFdens[:,0,:,:]
	1619	var1 = ncobj.variables[depvars[1]][:]
	1620	var2 = ncobj.variables[depvars[2]][:]
	1621	var3 = ncobj.variables[depvars[3]][:]
	1622	var4 = ncobj.variables[depvars[4]][:]
	1623	var5 = ncobj.variables[depvars[5]][:]
	1624	var6 = ncobj.variables[depvars[6]][:,0,:,:]
	1625
	1626	dnamesvar = list(ncobj.variables[depvars[1]].dimensions)
	1627	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1628
	1629	diagout = np.zeros(var1.shape, dtype=np.float)
	1630	diagout, diagoutd, diagoutvd = diag.Forcompute_potevap_orPM(var0, var1, var2,\
	1631	var3, var4, var5, var6, dnamesvar, dvnamesvar)
	1632
	1633	# Removing the nonChecking variable-dimensions from the initial list
	1634	varsadd = []
	1635	for nonvd in NONchkvardims:
	1636	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1637	varsadd.append(nonvd)
	1638
	1639	ncvar.insert_variable(ncobj, 'evspsblpot', diagout, diagoutd, diagoutvd, newnc)
	1640
[1795]	1641	# WRFmslp_ptarget sea-level pressure following ECMWF method as PSFC, HGT, WRFt, WRFp, ZNU, ZNW
	1642	elif diagn == 'WRFpsl_ecmwf':
	1643	var0 = ncobj.variables[depvars[0]][:]
	1644	var1 = ncobj.variables[depvars[1]][0,:,:]
	1645	var2 = WRFt[:,0,:,:]
	1646	var4 = WRFp[:,0,:,:]
	1647	var5 = ncobj.variables[depvars[4]][0,:]
	1648	var6 = ncobj.variables[depvars[5]][0,:]
	1649
	1650	# This is quite too appriximate!! passing pressure at half-levels to 2nd full
	1651	# level, using eta values at full (ZNW) and half (ZNU) mass levels
	1652	var3 = WRFp[:,0,:,:] + (var6[1] - var5[0])*(WRFp[:,1,:,:] - WRFp[:,0,:,:])/ \
	1653	(var5[1]-var5[0])
	1654
	1655	dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
	1656	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1657
	1658	diagout = np.zeros(var0.shape, dtype=np.float)
	1659	diagout, diagoutd, diagoutvd = diag.Forcompute_psl_ecmwf(var0, var1, var2, \
	1660	var3, var4, dnamesvar, dvnamesvar)
	1661
	1662	# Removing the nonChecking variable-dimensions from the initial list
	1663	varsadd = []
	1664	for nonvd in NONchkvardims:
	1665	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1666	varsadd.append(nonvd)
	1667
	1668	ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
	1669
[1758]	1670	# WRFmslp_ptarget sea-level pressure following ptarget method as WRFp, PSFC, WRFt, HGT, QVAPOR
	1671	elif diagn == 'WRFpsl_ptarget':
	1672	var0 = WRFp
	1673	var1 = ncobj.variables[depvars[1]][:]
	1674	var2 = WRFt
	1675	var3 = ncobj.variables[depvars[3]][0,:,:]
	1676	var4 = ncobj.variables[depvars[4]][:]
	1677
	1678	dnamesvar = list(ncobj.variables[depvars[4]].dimensions)
	1679	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1680
	1681	diagout = np.zeros(var0.shape, dtype=np.float)
	1682	diagout, diagoutd, diagoutvd = diag.Forcompute_psl_ptarget(var0, var1, var2, \
	1683	var3, var4, 700000., dnamesvar, dvnamesvar)
	1684
	1685	# Removing the nonChecking variable-dimensions from the initial list
	1686	varsadd = []
	1687	for nonvd in NONchkvardims:
	1688	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1689	varsadd.append(nonvd)
	1690
	1691	ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
	1692
[390]	1693	# WRFp pressure from WRF as P + PB
[1758]	1694	elif diagn == 'WRFp':
[1944]	1695	var0 = ncobj.variables[depvars[0]][:]
	1696	var1 = ncobj.variables[depvars[1]][:]
[365]	1697
[1944]	1698	diagout = var0 + var1
[1999]	1699	diagoutd = list(ncobj.variables[depvars[0]].dimensions)
	1700	diagoutvd = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
[365]	1701
[1999]	1702	# Removing the nonChecking variable-dimensions from the initial list
	1703	varsadd = []
	1704	for nonvd in NONchkvardims:
	1705	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1706	varsadd.append(nonvd)
[365]	1707
[1999]	1708	ncvar.insert_variable(ncobj, 'pres', diagout, diagoutd, diagoutvd, newnc)
	1709
[365]	1710	# WRFpos
[1758]	1711	elif diagn == 'WRFpos':
[365]	1712
	1713	dnamesvar = ncobj.variables['MAPFAC_M'].dimensions
	1714	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1715
	1716	ncvar.insert_variable(ncobj, 'WRFpos', WRFpos, dnamesvar, dvnamesvar, newnc)
	1717
	1718	# WRFprw WRF water vapour path WRFdens, QVAPOR
[1758]	1719	elif diagn == 'WRFprw':
[365]	1720
	1721	var0 = WRFdens
	1722	var1 = ncobj.variables[depvars[1]]
	1723
	1724	dnamesvar = list(var1.dimensions)
	1725	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1726
[1675]	1727	diagout, diagoutd, diagoutvd = diag.compute_prw(var0, var1, dnamesvar,dvnamesvar)
[365]	1728
[1586]	1729	# Removing the nonChecking variable-dimensions from the initial list
	1730	varsadd = []
	1731	diagoutvd = list(dvnames)
	1732	for nonvd in NONchkvardims:
	1733	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1734	varsadd.append(nonvd)
[365]	1735	ncvar.insert_variable(ncobj, 'prw', diagout, diagoutd, diagoutvd, newnc)
	1736
	1737	# WRFrh (P, T, QVAPOR)
[1758]	1738	elif diagn == 'WRFrh':
[365]	1739
	1740	dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
	1741	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1742
[878]	1743	ncvar.insert_variable(ncobj, 'hur', WRFrh, dnames, dvnames, newnc)
[365]	1744
	1745	# WRFrhs (PSFC, T2, Q2)
[1758]	1746	elif diagn == 'WRFrhs':
[365]	1747
	1748	var0 = ncobj.variables[depvars[0]][:]
	1749	var1 = ncobj.variables[depvars[1]][:]
	1750	var2 = ncobj.variables[depvars[2]][:]
	1751
	1752	dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
	1753	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1754
[1675]	1755	diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
[878]	1756	ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
[365]	1757
	1758	# rvors (u10, v10, WRFpos)
[1758]	1759	elif diagn == 'WRFrvors':
[365]	1760
	1761	var0 = ncobj.variables[depvars[0]]
	1762	var1 = ncobj.variables[depvars[1]]
	1763
	1764	diagout = rotational_z(var0, var1, distx)
	1765
	1766	dnamesvar = ncobj.variables[depvars[0]].dimensions
	1767	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1768
	1769	ncvar.insert_variable(ncobj, 'rvors', diagout, dnamesvar, dvnamesvar, newnc)
	1770
[884]	1771	# WRFt (T, P, PB)
[1758]	1772	elif diagn == 'WRFt':
[884]	1773	var0 = ncobj.variables[depvars[0]][:]
	1774	var1 = ncobj.variables[depvars[1]][:]
	1775	var2 = ncobj.variables[depvars[2]][:]
[654]	1776
[884]	1777	p0=100000.
	1778	p=var1 + var2
	1779
	1780	WRFt = (var0 + 300.)(p/p0)*(2./7.)
	1781
	1782	dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
	1783	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1784
[1382]	1785	# Removing the nonChecking variable-dimensions from the initial list
	1786	varsadd = []
[1389]	1787	diagoutvd = list(dvnames)
[1382]	1788	for nonvd in NONchkvardims:
[1389]	1789	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
[1382]	1790	varsadd.append(nonvd)
	1791
[1389]	1792	ncvar.insert_variable(ncobj, 'ta', WRFt, dnames, diagoutvd, newnc)
[884]	1793
[1942]	1794	# WRFtda (WRFrh, WRFt)
	1795	elif diagn == 'WRFtda':
	1796	ARM2 = fdef.module_definitions.arm2
	1797	ARM3 = fdef.module_definitions.arm3
	1798
	1799	gammatarh = np.log(WRFrh) + ARM2*(WRFt-273.15)/((WRFt-273.15)+ARM3)
[1943]	1800	td = ARM3*gammatarh/(ARM2-gammatarh)
[1942]	1801
	1802	dnamesvar = list(ncobj.variables['T'].dimensions)
	1803	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1804
	1805	# Removing the nonChecking variable-dimensions from the initial list
	1806	varsadd = []
	1807	diagoutvd = list(dvnames)
	1808	for nonvd in NONchkvardims:
	1809	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1810	varsadd.append(nonvd)
	1811
	1812	ncvar.insert_variable(ncobj, 'tda', td, dnames, diagoutvd, newnc)
	1813
[1966]	1814	# WRFtdas (PSFC, T2, Q2)
	1815	elif diagn == 'WRFtdas':
[1962]	1816	ARM2 = fdef.module_definitions.arm2
	1817	ARM3 = fdef.module_definitions.arm3
	1818
	1819	var0 = ncobj.variables[depvars[0]][:]
	1820	var1 = ncobj.variables[depvars[1]][:]
	1821	var2 = ncobj.variables[depvars[2]][:]
	1822
	1823	dnamesvar = list(ncobj.variables[depvars[1]].dimensions)
	1824	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1825
	1826	rhs, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
	1827
	1828	gammatarhs = np.log(rhs) + ARM2*(var1-273.15)/((var1-273.15)+ARM3)
[1970]	1829	tdas = ARM3*gammatarhs/(ARM2-gammatarhs) + 273.15
[1962]	1830
	1831	# Removing the nonChecking variable-dimensions from the initial list
	1832	varsadd = []
	1833	diagoutvd = list(dvnames)
	1834	for nonvd in NONchkvardims:
	1835	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1836	varsadd.append(nonvd)
	1837
[1966]	1838	ncvar.insert_variable(ncobj, 'tdas', tdas, dnames, diagoutvd, newnc)
[1962]	1839
[914]	1840	# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
[1758]	1841	elif diagn == 'WRFua':
[914]	1842	var0 = ncobj.variables[depvars[0]][:]
	1843	var1 = ncobj.variables[depvars[1]][:]
	1844	var2 = ncobj.variables[depvars[2]][:]
	1845	var3 = ncobj.variables[depvars[3]][:]
	1846
	1847	# un-staggering variables
[1999]	1848	if len(var0.shape) == 4:
	1849	unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
	1850	elif len(var0.shape) == 3:
	1851	# Asuming sunding point (dimt, dimz, dimstgx)
	1852	unstgdims = [var0.shape[0], var0.shape[1]]
	1853
[914]	1854	ua = np.zeros(tuple(unstgdims), dtype=np.float)
	1855	unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
	1856	unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
	1857
[1999]	1858	if len(var0.shape) == 4:
	1859	unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
	1860	unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
[914]	1861
[1999]	1862	for iz in range(var0.shape[1]):
	1863	ua[:,iz,:,:] = unstgvar0[:,iz,:,:]var3 - unstgvar1[:,iz,:,:]var2
	1864
	1865	dnamesvar = ['Time','bottom_top','south_north','west_east']
	1866
	1867	elif len(var0.shape) == 3:
	1868	unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1])
	1869	unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1])
	1870	for iz in range(var0.shape[1]):
	1871	ua[:,iz] = unstgvar0[:,iz]var3 - unstgvar1[:,iz]var2
	1872
	1873	dnamesvar = ['Time','bottom_top']
	1874
[914]	1875	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1876
[1404]	1877	# Removing the nonChecking variable-dimensions from the initial list
	1878	varsadd = []
	1879	diagoutvd = list(dvnames)
	1880	for nonvd in NONchkvardims:
	1881	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1882	varsadd.append(nonvd)
[914]	1883
[1404]	1884	ncvar.insert_variable(ncobj, 'ua', ua, dnames, diagoutvd, newnc)
	1885
[914]	1886	# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
[1758]	1887	elif diagn == 'WRFva':
[914]	1888	var0 = ncobj.variables[depvars[0]][:]
	1889	var1 = ncobj.variables[depvars[1]][:]
	1890	var2 = ncobj.variables[depvars[2]][:]
	1891	var3 = ncobj.variables[depvars[3]][:]
	1892
	1893	# un-staggering variables
[1999]	1894	if len(var0.shape) == 4:
	1895	unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
	1896	elif len(var0.shape) == 3:
	1897	# Asuming sunding point (dimt, dimz, dimstgx)
	1898	unstgdims = [var0.shape[0], var0.shape[1]]
	1899
[914]	1900	va = np.zeros(tuple(unstgdims), dtype=np.float)
	1901	unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
	1902	unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
	1903
[1999]	1904	if len(var0.shape) == 4:
	1905	unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
	1906	unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
	1907
	1908	for iz in range(var0.shape[1]):
	1909	va[:,iz,:,:] = unstgvar0[:,iz,:,:]var2 + unstgvar1[:,iz,:,:]var3
	1910
	1911	dnamesvar = ['Time','bottom_top','south_north','west_east']
	1912
	1913	elif len(var0.shape) == 3:
	1914	unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1])
	1915	unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1])
	1916	for iz in range(var0.shape[1]):
	1917	va[:,iz] = unstgvar0[:,iz]var2 + unstgvar1[:,iz]var3
	1918
	1919	dnamesvar = ['Time','bottom_top']
	1920
[914]	1921	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1922
[1404]	1923	# Removing the nonChecking variable-dimensions from the initial list
	1924	varsadd = []
	1925	diagoutvd = list(dvnames)
	1926	for nonvd in NONchkvardims:
	1927	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1928	varsadd.append(nonvd)
	1929	ncvar.insert_variable(ncobj, 'va', va, dnames, diagoutvd, newnc)
[914]	1930
[1980]	1931
	1932	# WRFwd (U, V, SINALPHA, COSALPHA) to be rotated !!
	1933	elif diagn == 'WRFwd':
	1934	var0 = ncobj.variables[depvars[0]][:]
	1935	var1 = ncobj.variables[depvars[1]][:]
	1936	var2 = ncobj.variables[depvars[2]][:]
	1937	var3 = ncobj.variables[depvars[3]][:]
	1938
	1939	# un-staggering variables
[1999]	1940	if len(var0.shape) == 4:
	1941	unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
	1942	elif len(var0.shape) == 3:
	1943	# Asuming sunding point (dimt, dimz, dimstgx)
	1944	unstgdims = [var0.shape[0], var0.shape[1]]
	1945
[1980]	1946	ua = np.zeros(tuple(unstgdims), dtype=np.float)
	1947	va = np.zeros(tuple(unstgdims), dtype=np.float)
	1948	unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
	1949	unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
	1950
[1999]	1951	if len(var0.shape) == 4:
	1952	unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
	1953	unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
[1980]	1954
[1999]	1955	for iz in range(var0.shape[1]):
	1956	ua[:,iz,:,:] = unstgvar0[:,iz,:,:]var3 - unstgvar1[:,iz,:,:]var2
	1957	va[:,iz,:,:] = unstgvar0[:,iz,:,:]var2 + unstgvar1[:,iz,:,:]var3
	1958
	1959	dnamesvar = ['Time','bottom_top','south_north','west_east']
	1960
	1961	elif len(var0.shape) == 3:
	1962	unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1])
	1963	unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1])
	1964	for iz in range(var0.shape[1]):
	1965	ua[:,iz] = unstgvar0[:,iz]var3 - unstgvar1[:,iz]var2
	1966	va[:,iz] = unstgvar0[:,iz]var2 + unstgvar1[:,iz]var3
	1967
	1968	dnamesvar = ['Time','bottom_top']
	1969
[1980]	1970	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	1971
	1972	wd, dnames, dvnames = diag.compute_wd(ua, va, dnamesvar, dvnamesvar)
	1973
	1974	# Removing the nonChecking variable-dimensions from the initial list
	1975	varsadd = []
	1976	diagoutvd = list(dvnames)
	1977	for nonvd in NONchkvardims:
	1978	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	1979	varsadd.append(nonvd)
	1980
	1981	ncvar.insert_variable(ncobj, 'wd', wd, dnames, diagoutvd, newnc)
	1982
[914]	1983	# WRFtime
[1758]	1984	elif diagn == 'WRFtime':
[654]	1985
	1986	diagout = WRFtime
	1987
	1988	dnamesvar = ['Time']
	1989	dvnamesvar = ['Times']
	1990
	1991	ncvar.insert_variable(ncobj, 'time', diagout, dnamesvar, dvnamesvar, newnc)
	1992
[959]	1993	# ws (U, V)
[1758]	1994	elif diagn == 'ws':
[959]	1995
	1996	# un-staggering variables
[1999]	1997	if len(var0.shape) == 4:
	1998	unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
	1999	elif len(var0.shape) == 3:
	2000	# Asuming sunding point (dimt, dimz, dimstgx)
	2001	unstgdims = [var0.shape[0], var0.shape[1]]
	2002
	2003	ua = np.zeros(tuple(unstgdims), dtype=np.float)
[959]	2004	va = np.zeros(tuple(unstgdims), dtype=np.float)
	2005	unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
	2006	unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
	2007
[1999]	2008	if len(var0.shape) == 4:
	2009	unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
	2010	unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
	2011
	2012	for iz in range(var0.shape[1]):
	2013	ua[:,iz,:,:] = unstgvar0[:,iz,:,:]var3 - unstgvar1[:,iz,:,:]var2
	2014	va[:,iz,:,:] = unstgvar0[:,iz,:,:]var2 + unstgvar1[:,iz,:,:]var3
	2015
	2016	dnamesvar = ['Time','bottom_top','south_north','west_east']
	2017
	2018	elif len(var0.shape) == 3:
	2019	unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1])
	2020	unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1])
	2021	for iz in range(var0.shape[1]):
	2022	ua[:,iz] = unstgvar0[:,iz]var3 - unstgvar1[:,iz]var2
	2023	va[:,iz] = unstgvar0[:,iz]var2 + unstgvar1[:,iz]var3
	2024
	2025	dnamesvar = ['Time','bottom_top']
	2026
[959]	2027	diagout = np.sqrt(unstgvar0unstgvar0 + unstgvar1unstgvar1)
	2028
	2029	# dnamesvar = ncobj.variables[depvars[0]].dimensions
	2030	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	2031
[1408]	2032	# Removing the nonChecking variable-dimensions from the initial list
	2033	varsadd = []
	2034	diagoutvd = list(dvnamesvar)
	2035	for nonvd in NONchkvardims:
	2036	if gen.searchInlist(dvnamesvar,nonvd): diagoutvd.remove(nonvd)
	2037	varsadd.append(nonvd)
	2038	ncvar.insert_variable(ncobj, 'ws', diagout, dnamesvar, diagoutvd, newnc)
[959]	2039
[365]	2040	# wss (u10, v10)
[1758]	2041	elif diagn == 'wss':
[365]	2042
	2043	var0 = ncobj.variables[depvars[0]][:]
	2044	var1 = ncobj.variables[depvars[1]][:]
	2045
	2046	diagout = np.sqrt(var0var0 + var1var1)
	2047
	2048	dnamesvar = ncobj.variables[depvars[0]].dimensions
	2049	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
	2050
	2051	ncvar.insert_variable(ncobj, 'wss', diagout, dnamesvar, dvnamesvar, newnc)
	2052
[654]	2053	# WRFheight height from WRF geopotential as WRFGeop/g
[1758]	2054	elif diagn == 'WRFheight':
[654]	2055
	2056	diagout = WRFgeop/grav
	2057
[1412]	2058	# Removing the nonChecking variable-dimensions from the initial list
	2059	varsadd = []
	2060	diagoutvd = list(dvnames)
	2061	for nonvd in NONchkvardims:
	2062	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	2063	varsadd.append(nonvd)
[654]	2064
[1412]	2065	ncvar.insert_variable(ncobj, 'zhgt', diagout, dnames, diagoutvd, newnc)
	2066
[1413]	2067	# WRFheightrel relative-height from WRF geopotential as WRFgeop(PH + PHB)/g-HGT 'WRFheightrel\|PH@PHB@HGT
[1758]	2068	elif diagn == 'WRFheightrel':
[1413]	2069	var0 = ncobj.variables[depvars[0]][:]
	2070	var1 = ncobj.variables[depvars[1]][:]
	2071	var2 = ncobj.variables[depvars[2]][:]
	2072
	2073	dimz = var0.shape[1]
	2074	diagout = np.zeros(tuple(var0.shape), dtype=np.float)
	2075	for iz in range(dimz):
[1419]	2076	diagout[:,iz,:,:] = (var0[:,iz,:,:]+ var1[:,iz,:,:])/grav - var2
[1413]	2077
	2078	# Removing the nonChecking variable-dimensions from the initial list
	2079	varsadd = []
	2080	diagoutvd = list(dvnames)
	2081	for nonvd in NONchkvardims:
	2082	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	2083	varsadd.append(nonvd)
	2084
	2085	ncvar.insert_variable(ncobj, 'zhgtrel', diagout, dnames, diagoutvd, newnc)
	2086
[1773]	2087	# WRFzmla_gen generic boundary layer hieght computation from WRF theta, QVAPOR, WRFgeop, HGT,
	2088	elif diagn == 'WRFzmlagen':
	2089	var0 = ncobj.variables[depvars[0]][:]+300.
	2090	var1 = ncobj.variables[depvars[1]][:]
	2091	dimz = var0.shape[1]
	2092	var2 = WRFgeop[:,1:dimz+1,:,:]/9.8
	2093	var3 = ncobj.variables[depvars[3]][0,:,:]
	2094
	2095	diagout, diagoutd, diagoutvd = diag.Forcompute_zmla_gen(var0,var1,var2,var3, \
	2096	dnames,dvnames)
	2097
	2098	# Removing the nonChecking variable-dimensions from the initial list
	2099	varsadd = []
	2100	for nonvd in NONchkvardims:
	2101	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	2102	varsadd.append(nonvd)
	2103
	2104	ncvar.insert_variable(ncobj, 'zmla', diagout, diagoutd, diagoutvd, newnc)
	2105
[1784]	2106	# WRFzwind wind extrapolation at a given height using power law computation from WRF
	2107	# U, V, WRFz, U10, V10, SINALPHA, COSALPHA, z=[zval]
[1776]	2108	elif diagn == 'WRFzwind':
	2109	var0 = ncobj.variables[depvars[0]][:]
	2110	var1 = ncobj.variables[depvars[1]][:]
[1777]	2111	var2 = WRFz
[1776]	2112	var3 = ncobj.variables[depvars[3]][:]
	2113	var4 = ncobj.variables[depvars[4]][:]
	2114	var5 = ncobj.variables[depvars[5]][0,:,:]
	2115	var6 = ncobj.variables[depvars[6]][0,:,:]
[1777]	2116	var7 = np.float(depvars[7].split('=')[1])
[1776]	2117
	2118	# un-staggering 3D winds
	2119	unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
	2120	va = np.zeros(tuple(unstgdims), dtype=np.float)
	2121	unvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
	2122	unvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
	2123	unvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
	2124	unvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
	2125
	2126	diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwind(unvar0, \
[1777]	2127	unvar1, var2, var3, var4, var5, var6, var7, dnames, dvnames)
[1776]	2128
	2129	# Removing the nonChecking variable-dimensions from the initial list
	2130	varsadd = []
	2131	for nonvd in NONchkvardims:
	2132	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	2133	varsadd.append(nonvd)
	2134
	2135	ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
	2136	ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)
	2137
[1784]	2138	# WRFzwind wind extrapolation at a given hieght using logarithmic law computation
	2139	# from WRF U, V, WRFz, U10, V10, SINALPHA, COSALPHA, z=[zval]
	2140	elif diagn == 'WRFzwind_log':
	2141	var0 = ncobj.variables[depvars[0]][:]
	2142	var1 = ncobj.variables[depvars[1]][:]
	2143	var2 = WRFz
	2144	var3 = ncobj.variables[depvars[3]][:]
	2145	var4 = ncobj.variables[depvars[4]][:]
	2146	var5 = ncobj.variables[depvars[5]][0,:,:]
	2147	var6 = ncobj.variables[depvars[6]][0,:,:]
	2148	var7 = np.float(depvars[7].split('=')[1])
	2149
	2150	# un-staggering 3D winds
	2151	unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
	2152	va = np.zeros(tuple(unstgdims), dtype=np.float)
	2153	unvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
	2154	unvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
	2155	unvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
	2156	unvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
	2157
	2158	diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwind_log(unvar0, \
	2159	unvar1, var2, var3, var4, var5, var6, var7, dnames, dvnames)
	2160
	2161	# Removing the nonChecking variable-dimensions from the initial list
	2162	varsadd = []
	2163	for nonvd in NONchkvardims:
	2164	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	2165	varsadd.append(nonvd)
	2166
	2167	ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
	2168	ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)
	2169
[1783]	2170	# WRFzwindMO wind extrapolation at a given height computation using Monin-Obukhow
	2171	# theory from WRF UST, ZNT, RMOL, U10, V10, SINALPHA, COSALPHA, z=[zval]
[1784]	2172	# NOTE: only useful for [zval] < 80. m
[1783]	2173	elif diagn == 'WRFzwindMO':
	2174	var0 = ncobj.variables[depvars[0]][:]
	2175	var1 = ncobj.variables[depvars[1]][:]
	2176	var2 = ncobj.variables[depvars[2]][:]
	2177	var3 = ncobj.variables[depvars[3]][:]
	2178	var4 = ncobj.variables[depvars[4]][:]
	2179	var5 = ncobj.variables[depvars[5]][0,:,:]
	2180	var6 = ncobj.variables[depvars[6]][0,:,:]
	2181	var7 = np.float(depvars[7].split('=')[1])
	2182
	2183	diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwindMO(var0, var1,\
	2184	var2, var3, var4, var5, var6, var7, dnames, dvnames)
	2185
	2186	# Removing the nonChecking variable-dimensions from the initial list
	2187	varsadd = []
	2188	for nonvd in NONchkvardims:
	2189	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	2190	varsadd.append(nonvd)
	2191
	2192	ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
	2193	ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)
	2194
[2619]	2195	# zmla_gen generic boundary layer hieght computation from generic 2D-file theta, qv, zg, orog,
	2196	elif diagn == 'zmlagen':
	2197	var0 = ncobj.variables[depvars[0]][:]
	2198	var1 = ncobj.variables[depvars[1]][:]
	2199	dimz = var0.shape[1]
	2200	var2 = ncobj.variables[depvars[2]][:]
	2201	var3 = ncobj.variables[depvars[3]][:]
	2202
	2203	diagout, diagoutd, diagoutvd = diag.Forcompute_zmla_gen(var0,var1,var2,var3, \
	2204	dnames,dvnames)
	2205
	2206	# Removing the nonChecking variable-dimensions from the initial list
	2207	varsadd = []
	2208	for nonvd in NONchkvardims:
	2209	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	2210	varsadd.append(nonvd)
	2211
	2212	ncvar.insert_variable(ncobj, 'zmla', diagout, diagoutd, diagoutvd, newnc)
	2213
	2214	# zmla_genUWsnd generic boundary layer hieght computation from UWyoming sounding file theta, qv, zg
	2215	elif diagn == 'zmlagenUWsnd':
	2216	var0 = ncobj.variables[depvars[0]][:]
	2217	var1 = ncobj.variables[depvars[1]][:]
	2218	dimz = var0.shape[1]
[2620]	2219	var2 = ncobj.variables[depvars[2]][:]/9.8
[2619]	2220	var3 = ncobj.getncattr('Station_elevation')
	2221
	2222	diagout, diagoutd, diagoutvd = diag.Forcompute_zmla_gen(var0,var1,var2,var3, \
	2223	dnames,dvnames)
	2224
[2620]	2225	# No need to remove topography height
	2226	diagout = diagout + var3
	2227
[2619]	2228	# Removing the nonChecking variable-dimensions from the initial list
	2229	varsadd = []
	2230	for nonvd in NONchkvardims:
	2231	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
	2232	varsadd.append(nonvd)
	2233
	2234	ncvar.insert_variable(ncobj, 'zmla', diagout, diagoutd, diagoutvd, newnc)
	2235
[365]	2236	else:
	2237	print errormsg
[1758]	2238	print ' ' + main + ": diagnostic '" + diagn + "' not ready!!!"
[365]	2239	print ' available diagnostics: ', availdiags
	2240	quit(-1)
	2241
	2242	newnc.sync()
[1351]	2243	# Adding that additional variables required to compute some diagnostics which
	2244	# where not in the original file
[1944]	2245	print ' adding additional variables...'
[1351]	2246	for vadd in varsadd:
[1942]	2247	if not gen.searchInlist(newnc.variables.keys(),vadd) and \
	2248	dictcompvars.has_key(vadd):
[1351]	2249	attrs = dictcompvars[vadd]
	2250	vvn = attrs['name']
	2251	if not gen.searchInlist(newnc.variables.keys(), vvn):
	2252	iidvn = dvnames.index(vadd)
	2253	dnn = dnames[iidvn]
	2254	if vadd == 'WRFtime':
	2255	dvarvals = WRFtime[:]
	2256	newvar = newnc.createVariable(vvn, 'f8', (dnn))
	2257	newvar[:] = dvarvals
	2258	for attn in attrs.keys():
	2259	if attn != 'name':
	2260	attv = attrs[attn]
	2261	ncvar.set_attribute(newvar, attn, attv)
[365]	2262
	2263	# end of diagnostics
	2264
	2265	# Global attributes
	2266	##
[1758]	2267	ncvar.add_global_PyNCplot(newnc, main, None, '2.0')
[365]	2268
	2269	gorigattrs = ncobj.ncattrs()
	2270	for attr in gorigattrs:
	2271	attrv = ncobj.getncattr(attr)
	2272	atvar = ncvar.set_attribute(newnc, attr, attrv)
	2273
	2274	ncobj.close()
	2275	newnc.close()
	2276
	2277	print '\n' + main + ': successfull writting of diagnostics file "' + ofile + '" !!!'

Note: See TracBrowser for help on using the repository browser.

Download in other formats: