Context Navigation

diagnostics.py @ 2788

Last change on this file since 2788 was 2767, checked in by lfita, 6 years ago

Adding:

`gradient2Dh': 1st order Horizontal 2D-gradient of any variable

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

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

Context Navigation

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

Download in other formats: