Context Navigation

diagnostics.py @ 2828

Last change on this file since 2828 was 2767, checked in by lfita, 5 years ago

Adding:

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

File size: 106.3 KB

Line
1	# Python script to comput diagnostics
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
10	# File diagnostics.inf provides the combination of variables to get the desired diagnostic
11	# To be used with module_ForDiagnostics.F90, module_ForDiagnosticsVars.F90, module_generic.F90
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
14
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
16	## e.g. # diagnostics.py -f /home/lluis/PY/diagnostics.inf -d variable_combo -v WRFprc
17
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)
29	# compute_range_faces: Function to compute faces [uphill, valley, downhill] of sections of a mountain
30	# range, along a given face
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
41	# Forcompute_psl_ptarget: Function to compute the sea-level pressure following target_pressure value found in `p_interp.F'
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
51	# timeoverthres: When a given variable [varname] overpass a given [value]. Being [CFvarn] the name of the diagnostics in
52	# variables_values.dat
53	# timemax ([varname], time). When a given variable [varname] got its maximum
54
55	from optparse import OptionParser
56	import numpy as np
57	import numpy.ma as ma
58	from netCDF4 import Dataset as NetCDFFile
59	import os
60	import re
61	import nc_var_tools as ncvar
62	import generic_tools as gen
63	import datetime as dtime
64	import module_ForDiag as fdin
65	import module_ForDef as fdef
66	import diag_tools as diag
67
68	main = 'diagnostics.py'
69	errormsg = 'ERROR -- error -- ERROR -- error'
70	warnmsg = 'WARNING -- warning -- WARNING -- warning'
71
72	# Constants
73	grav = 9.81
74
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",
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,
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	#######
92	availdiags = ['ACRAINTOT', 'accum', 'clt', 'cllmh', 'convini', 'deaccum', 'fog_K84', \
93	'fog_RUC', 'front_R04', 'frontogenesis', 'gradient2Dh', 'rhs_tas_tds', 'LMDZrh', \
94	'mslp', 'OMEGAw', 'RAINTOT', 'range_faces', \
95	'rvors', 'td', 'timemax', 'timeoverthres', 'turbulence', 'tws', 'uavaFROMwswd', \
96	'WRFbnds', 'WRFcape_afwa', 'WRFclivi', 'WRFclwvi', 'WRF_denszint', 'WRFgeop', \
97	'WRFmrso', 'WRFmrsos', 'WRFpotevap_orPM', 'WRFp', 'WRFpsl_ecmwf', \
98	'WRFpsl_ptarget', 'WRFrvors', 'WRFslw', 'ws', 'wds', 'wss', 'WRFheight', \
99	'WRFheightrel', 'WRFtda', 'WRFtdas', 'WRFtws', 'WRFua', 'WRFva', 'WRFzwind', \
100	'WRFzwind_log', 'WRFzwindMO', 'zmlagen', 'zmlagenUWsnd']
101
102	methods = ['accum', 'deaccum']
103
104	# Variables not to check
105	NONcheckingvars = ['accum', 'cllmh', 'deaccum', 'face', 'LONLATdxdy', 'LONLATdx', \
106	'LONLATdy', 'params', 'reglonlatbnds', 'TSrhs', 'TStd', 'TSwds', 'TSwss', \
107	'UNua', 'UNva', 'UNwa', \
108	'WRFbils', 'WRFbnds', \
109	'WRFclivi', 'WRFclwvi', 'WRFdens', 'WRFdx', 'WRFdxdy', 'WRFdxdywps', 'WRFdy', \
110	'WRFdz', 'WRFgeop', 'WRFp', 'WRFtd', \
111	'WRFpos', 'WRFprc', 'WRFprls', 'WRFrh', 'LMDZrh', 'LMDZrhs', \
112	'WRFrhs', 'WRFrvors', \
113	'WRFt', 'WRFtime', 'WRFua', 'WRFva', 'WRFwds', 'WRFwss', 'WRFheight', 'WRFz', \
114	'WRFzg']
115
116	# diagnostics not to check their dependeny
117	NONcheckdepvars = ['accum', 'deaccum', 'timeoverthres', 'WRF_denszint', \
118	'WRFzwind_log', 'WRFzwind', 'WRFzwindMO']
119
120	NONchkvardims = ['WRFtime']
121
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
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
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
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
170	WRFz_compute = False
171	WRFdx_compute = False
172	WRFdy_compute = False
173	WRFdz_compute = False
174	WRFdxdy_compute = False
175	WRFdxdywps_compute = False
176	LONLATdxdy_compute = False
177	LONLATdx_compute = False
178	LONLATdy_compute = False
179	UNua_compute = False
180	UNva_compute = False
181	UNwa_compute = False
182
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:
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
205	if dnv == 'WRFz':WRFz_compute = True
206	if dnv == 'WRFdx':WRFdx_compute = True
207	if dnv == 'WRFdy':WRFdy_compute = True
208	if dnv == 'WRFdz':WRFdz_compute = True
209	if dnv == 'WRFdxdy':WRFdxdy_compute = True
210	if dnv == 'WRFdxdywps':WRFdxdywps_compute = True
211	if dnv == 'LONLATdxdy':LONLATdxdy_compute = True
212	if dnv == 'LONLATdx':LONLATdx_compute = True
213	if dnv == 'LONLATdy':LONLATdy_compute = True
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]
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]
231	if depvars == 'WRFgeop':WRFgeop_compute = True
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
238	if depvars == 'WRFtime': WRFtime_compute = True
239	if depvars == 'WRFz': WRFz_compute = True
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
246	if depvars == 'LONLATdx': LONLATdx_compute = True
247	if depvars == 'LONLATdy': LONLATdy_compute = True
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]
257	else:
258	depvars = diags[idiag].split('\|')[1].split('@')
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
267	if gen.searchInlist(depvars, 'WRFz'): WRFz_compute = True
268	if gen.searchInlist(depvars, 'WRFdx'): WRFdx_compute = True
269	if gen.searchInlist(depvars, 'WRFdy'): WRFdy_compute = True
270	if gen.searchInlist(depvars, 'WRFdz'): WRFdz_compute = True
271	if gen.searchInlist(depvars, 'WRFdxdy'): WRFdxdy_compute = True
272	if gen.searchInlist(depvars, 'WRFdxdywps'): WRFdxdywps_compute = True
273	if gen.searchInlist(depvars, 'LONLATdxdy'): LONLATdxdy_compute = True
274	if gen.searchInlist(depvars, 'LONLATdx'): LONLATdx_compute = True
275	if gen.searchInlist(depvars, 'LONLATdy'): LONLATdy_compute = True
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]
291
292	# Dictionary with the new computed variables to be able to add them
293	dictcompvars = {}
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
299	# Attributes of the variable
300	Vvals = gen.variables_values('WRFgeop')
301	dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
302	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
303
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
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
314	if WRFght_compute:
315	print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...'
316	WRFght = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
317
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
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
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
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
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
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
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
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
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
431	if len(timeobj.shape) == 2:
432	dt = timeobj.shape[0]
433	else:
434	dt = 1
435	WRFtime = np.zeros((dt), dtype=np.float)
436
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')
445	WRFtime[0] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
446
447	tunits = tunitsval + ' since ' + refdateS
448
449	# Attributes of the variable
450	dictcompvars['WRFtime'] = {'name': 'time', 'standard_name': 'time', \
451	'long_name': 'time', 'units': tunits, 'calendar': 'gregorian'}
452
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')
469	dictcompvars['WRFz'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
470	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
471
472	if WRFdxdy_compute or WRFdx_compute or WRFdy_compute:
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')
501	dictcompvars['WRFds'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
502	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
503
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
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))'
539	PH = ncobj.variables['PH'][:]
540	PHB = ncobj.variables['PHB'][:]
541
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
559	if LONLATdxdy_compute or LONLATdx_compute or LONLATdy_compute :
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')
587	dictcompvars['WRFdx'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
588	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
589	Vvals = gen.variables_values('WRFdy')
590	dictcompvars['WRFdy'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
591	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
592	Vvals = gen.variables_values('WRFds')
593	dictcompvars['WRFds'] = {'name': Vvals[0], 'standard_name': Vvals[1], \
594	'long_name': Vvals[4].replace('\|',' '), 'units': Vvals[5]}
595
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
659	### ## #
660	# Going for the diagnostics
661	### ## #
662	print ' ' + main + ' ...'
663	varsadd = []
664
665	Varns = ncobj.variables.keys()
666	Varns.sort()
667
668	for idiag in range(Ndiags):
669	print ' diagnostic:',diags[idiag]
670	diagn = diags[idiag].split('\|')[0]
671	depvars = diags[idiag].split('\|')[1].split('@')
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:
678	# Checking without extra arguments of a depending variable (':', separated)
679	if depv.find(':') != -1: depv=depv.split(':')[0]
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 + "' !!"
687	print ' available ones:', Varns
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):
694	print errormsg
695	print ' ' + main + ": file '" + opts.ncfile + \
696	"' does not have variable '" + depvars + "' !!"
697	print ' available ones:', Varns
698	quit(-1)
699
700	print "\n Computing '" + diagn + "' from: ", depvars, '...'
701
702	# acraintot: accumulated total precipitation from WRF RAINC, RAINNC, RAINSH
703	if diagn == 'ACRAINTOT':
704
705	var0 = ncobj.variables[depvars[0]]
706	var1 = ncobj.variables[depvars[1]]
707	var2 = ncobj.variables[depvars[2]]
708
709	diagout = var0[:] + var1[:] + var2[:]
710
711	dnamesvar = var0.dimensions
712	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
713
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
720	ncvar.insert_variable(ncobj, 'pracc', diagout, dnamesvar, dvnamesvar, newnc)
721
722	# accum: acumulation of any variable as (Variable, time [as [tunits]
723	# from/since ....], newvarname)
724	elif diagn == 'accum':
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
732	diagout, diagoutd, diagoutvd = diag.compute_accum(var0,dnamesvar,dvnamesvar)
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)
739
740	CFvarn = gen.variables_values(depvars[0])[0]
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
750	dtime = (var1[1] - var1[0])*diag.timeunits_seconds(tunits)
751
752	ncvar.insert_variable(ncobj, CFvarn + 'acc', diagout*dtime, diagoutd, diagoutvd, newnc)
753
754	# cllmh with cldfra, pres
755	elif diagn == 'cllmh':
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
767	diagout, diagoutd, diagoutvd = diag.Forcompute_cllmh(var0,var1,dnames,dvnames)
768
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
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
780	elif diagn == 'clt':
781
782	var0 = ncobj.variables[depvars]
783	diagout, diagoutd, diagoutvd = diag.Forcompute_clt(var0,dnames,dvnames)
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
791	ncvar.insert_variable(ncobj, 'clt', diagout, diagoutd, diagoutvd, newnc)
792
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
814	# deaccum: deacumulation of any variable as (Variable, time [as [tunits]
815	# from/since ....], newvarname)
816	elif diagn == 'deaccum':
817
818	var0 = ncobj.variables[depvars[0]]
819	var1 = ncobj.variables[depvars[1]]
820
821	dnamesvar = var0.dimensions
822	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
823
824	diagout, diagoutd, diagoutvd = diag.compute_deaccum(var0,dnamesvar,dvnamesvar)
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)
831
832	# Transforming to a flux
833	if depvars[1] == 'XTIME':
834	dtimeunits = var1.getncattr('description')
835	tunits = dtimeunits.split(' ')[0]
836	else:
837	dtimeunits = var1.getncattr('units')
838	tunits = dtimeunits.split(' ')[0]
839
840	dtime = (var1[1] - var1[0])*diag.timeunits_seconds(tunits)
841	ncvar.insert_variable(ncobj, depvars[2], diagout/dtime, diagoutd, diagoutvd, \
842	newnc)
843
844	# fog_K84: Computation of fog and visibility following Kunkel, (1984) as QCLOUD, QICE
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
864	# fog_RUC: Computation of fog and visibility following Kunkel, (1984) as QVAPOR,
865	# WRFt, WRFp or Q2, T2, PSFC
866	elif diagn == 'fog_RUC':
867
868	var0 = ncobj.variables[depvars[0]]
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]]
878
879	dnamesvar = list(var0.dimensions)
880	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
881
882	diag1, diag2, diagoutd, diagoutvd = diag.Forcompute_fog_RUC(var0, var1, var2,\
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
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
922	# front_R04: Computation of the presence of a front as defined by Rodrigues et al.
923	# (2004) (tas, uas, vas, dxs, dys, 'params',[dtas],[dwss])
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]]
931	if depvars[3] == 'WRFdx' or depvars[3] == 'LONLATdx': var3 = WRFdx
932	else: var3 = ncobj.variables[depvars[3]]
933	if depvars[4] == 'WRFdy' or depvars[4] == 'LONLATdy': var4 = WRFdy
934	else: var4 = ncobj.variables[depvars[4]]
935	par1 = np.float(depvars[5].split(':')[1])
936	par2 = np.float(depvars[5].split(':')[2])
937
938	dnamesvar = list(var0.dimensions)
939	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
940
941	diag1, diag2, diag3, diag4, diagoutd, diagoutvd = diag.Forcompute_front_R04( \
942	var0, var1, var2, var3, var4, par1, par2, dnamesvar, dvnamesvar)
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)
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])
970	ncobj.sync()
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')
980	ncobj.sync()
981	ncvar.insert_variable(ncobj, 'dt2', diag4, diagoutd, diagoutvd, newnc, \
982	fill=gen.fillValueF)
983	ovar = newnc.variables['dt2']
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])
990	ncobj.sync()
991
992	# frontogenesis (theta, ua, va, wa, press, dxs, dys, dzs, time)
993	elif diagn == 'frontogenesis':
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]
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]][:]
1006	if depvars[4] == 'WRFp': var4 = WRFp
1007	else: var4 = ncobj.variables[depvars[4]][:]
1008	if depvars[5] == 'WRFdx': var5 = WRFdx
1009	else: var5 = ncobj.variables[depvars[5]][:]
1010	if depvars[6] == 'WRFdy': var6 = WRFdy
1011	else: var6 = ncobj.variables[depvars[6]][:]
1012	if depvars[7] == 'WRFdz': var7 = WRFdz[0,0:dz,0,0]
1013	else: var7 = ncobj.variables[depvars[7]][:]
1014	if depvars[8] == 'WRFtime': var08 = WRFtime
1015	else: var08 = ncobj.variables[depvars[8]][:]
1016
1017	# Assuming monotonic time-axis...
1018	var8 = var08[1] - var08[0]
1019
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
1029
1030	# Unifying...
1031	#var5 = var5/var5
1032	#var6 = var6/var6
1033	#var7 = var7/var7
1034
1035	diag1, diag2, diag3, diag4, diag5, diag6, diag7, diag8, diag9, diag10, \
1036	diagoutd, diagoutvd = diag.Forcompute_frontogenesis(var0, var1, var2, var3,\
1037	var4, var5, var6, var7, var8, dnames, dvnames)
1038
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
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
1058	ncvar.insert_variable(ncobj, 'diabh', diag1, diagoutd, diagoutvd, newnc, \
1059	gen.fillValueF)
1060	newnc.renameVariable('diabh', 'xdiabh')
1061	newnc.sync()
1062	ovar = newnc.variables['xdiabh']
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)
1072	newnc.renameVariable('diabh', 'ydiabh')
1073	newnc.sync()
1074	ovar = newnc.variables['ydiabh']
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)
1084	newnc.renameVariable('diabh', 'zdiabh')
1085	newnc.sync()
1086	ovar = newnc.variables['zdiabh']
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)
1097	newnc.renameVariable('def', 'xdef')
1098	newnc.sync()
1099	ovar = newnc.variables['xdef']
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)
1109	newnc.renameVariable('def', 'ydef')
1110	newnc.sync()
1111	ovar = newnc.variables['ydef']
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)
1121	newnc.renameVariable('def', 'zdef')
1122	newnc.sync()
1123	ovar = newnc.variables['zdef']
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)
1134	newnc.renameVariable('tilt', 'xtilt')
1135	newnc.sync()
1136	ovar = newnc.variables['xtilt']
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)
1147	newnc.renameVariable('tilt', 'ytilt')
1148	newnc.sync()
1149	ovar = newnc.variables['ytilt']
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)
1160	newnc.renameVariable('div', 'zdiv')
1161	newnc.sync()
1162	ovar = newnc.variables['zdiv']
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)
1174	ncvar.basicvardef(newvar, 'component', 'axis component: x, y, z', '-')
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
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]]
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,:,:]
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
1209	CFvarn = gen.variables_values(depvars[0])[0]
1210	ncvar.insert_variable(ncobj, depvars[0], diagout, diagoutd, diagoutvd, \
1211	newnc, gen.fillValueF)
1212	newnc.sync()
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'
1217	varu = ovar.units
1218	ovar.setncattr('units', varu+'ds-1')
1219	newnc.sync()
1220
1221	# LMDZrh (pres, t, r)
1222	elif diagn == 'LMDZrh':
1223
1224	var0 = ncobj.variables[depvars[0]][:]
1225	var1 = ncobj.variables[depvars[1]][:]
1226	var2 = ncobj.variables[depvars[2]][:]
1227
1228	diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnames,dvnames)
1229	ncvar.insert_variable(ncobj, 'hur', diagout, diagoutd, diagoutvd, newnc)
1230
1231	# LMDZrhs (psol, t2m, q2m)
1232	elif diagn == 'LMDZrhs':
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
1241	diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
1242
1243	ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
1244
1245	# range_faces: LON, LAT, HGT, WRFdxdy, 'face:['WE'/'SN']:[dsfilt]:[dsnewrange]:[hvalleyrange]'
1246	elif diagn == 'range_faces':
1247
1248	var0 = ncobj.variables[depvars[0]][:]
1249	var1 = ncobj.variables[depvars[1]][:]
1250	var2 = ncobj.variables[depvars[2]][:]
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])
1255
1256	dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
1257	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
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,:,:]
1263	dnamesvar.pop(0)
1264	dvnamesvar.pop(0)
1265	else:
1266	hgt = var2[:]
1267
1268	orogmax, ptorogmax, dhgt, peaks, valleys, origfaces, diagout, diagoutd, \
1269	diagoutvd, rng, rngorogmax, ptrngorogmax= diag.Forcompute_range_faces(lon, \
1270	lat, hgt, WRFdx, WRFdy, WRFds, face, dsfilt, dsnewrange, hvalleyrange, \
1271	dnamesvar, dvnamesvar)
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)
1279
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
1284	# adding variables to output file
1285	if face == 'WE': axis = 'lon'
1286	elif face == 'SN': axis = 'lat'
1287
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
1293	ncvar.insert_variable(ncobj, 'orogmax', rngorogmax, diagoutd, diagoutvd, \
1294	newnc, fill=gen.fillValueF)
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
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
1323	ncvar.insert_variable(ncobj, 'orogderiv', dhgt, diagoutd, diagoutvd, newnc)
1324	ovar = newnc.variables['orogderiv']
1325	ncvar.set_attribute(ovar, 'deriv', axis)
1326
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)
1335	ncvar.set_attributek(ovar, 'dist_filter', dsfilt, 'F')
1336
1337	ncvar.insert_variable(ncobj, 'rangefaces', origfaces, diagoutd, diagoutvd, \
1338	newnc, fill=gen.fillValueI)
1339	ovar = newnc.variables['rangefaces']
1340	ncvar.set_attribute(ovar, 'face', face)
1341	ncvar.set_attributek(ovar, 'dist_newrange', dsnewrange, 'F')
1342	ncvar.set_attributek(ovar, 'h_valley_newrange', hvalleyrange, 'F')
1343
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
1381	# tws: Bet Wulb temperature following Stull 2011 (tas, hurs)
1382	elif diagn == 'tws':
1383
1384	var0 = ncobj.variables[depvars[0]][:]
1385	var1 = ncobj.variables[depvars[1]][:]
1386
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)
1393	ncvar.insert_variable(ncobj, 'tws', diagout, diagoutd, diagoutvd, newnc)
1394
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 = []
1405	dnamesvar.append('bnds')
1406	dnamesvar.append(ncobj.variables[depvars[0]].dimensions[1])
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)
1423	newnc.sync()
1424	ncvar.insert_variable(ncobj, 'lat_bnds', cellbndsy, diagoutd, diagoutvd, newnc)
1425	newnc.sync()
1426
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)
1442
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
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
1495	# mslp: mean sea level pressure (pres, psfc, terrain, temp, qv)
1496	elif diagn == 'mslp' or diagn == 'WRFmslp':
1497
1498	var1 = ncobj.variables[depvars[1]][:]
1499	var2 = ncobj.variables[depvars[2]][:]
1500	var4 = ncobj.variables[depvars[4]][:]
1501
1502	if diagn == 'WRFmslp':
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
1513	diagout, diagoutd, diagoutvd = diag.compute_mslp(var0, var1, var2, var3, var4, \
1514	dnamesvar, dvnamesvar)
1515
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)
1522	ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
1523
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
1546	# OMEGAw (omega, p, t) from NCL formulation (https://www.ncl.ucar.edu/Document/Functions/Contributed/omega_to_w.shtml)
1547	elif diagn == 'OMEGAw':
1548
1549	var0 = ncobj.variables[depvars[0]][:]
1550	var1 = ncobj.variables[depvars[1]][:]
1551	var2 = ncobj.variables[depvars[2]][:]
1552
1553	dnamesvar = ncobj.variables[depvars[0]].dimensions
1554	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1555
1556	diagout, diagoutd, diagoutvd = diag.compute_OMEGAw(var0,var1,var2,dnamesvar,dvnamesvar)
1557
1558	ncvar.insert_variable(ncobj, 'wa', diagout, diagoutd, diagoutvd, newnc)
1559
1560	# raintot: instantaneous total precipitation from WRF as (RAINC + RAINC + RAINSH) / dTime
1561	elif diagn == 'RAINTOT':
1562
1563	var0 = ncobj.variables[depvars[0]]
1564	var1 = ncobj.variables[depvars[1]]
1565	var2 = ncobj.variables[depvars[2]]
1566
1567	if depvars[3] != 'WRFtime':
1568	var3 = ncobj.variables[depvars[3]]
1569	else:
1570	var3 = np.arange(var0.shape[0], dtype=int)
1571
1572	var = var0[:] + var1[:] + var2[:]
1573
1574	dnamesvar = var0.dimensions
1575	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1576
1577	diagout, diagoutd, diagoutvd = diag.compute_deaccum(var,dnamesvar,dvnamesvar)
1578
1579	# Transforming to a flux
1580	if var3.shape[0] > 1:
1581	if depvars[3] != 'WRFtime':
1582	dtimeunits = var3.getncattr('units')
1583	tunits = dtimeunits.split(' ')[0]
1584
1585	dtime = (var3[1] - var3[0])*diag.timeunits_seconds(tunits)
1586	else:
1587	var3 = ncobj.variables['Times']
1588	time1 = var3[0,:]
1589	time2 = var3[1,:]
1590	tmf1 = ''
1591	tmf2 = ''
1592	for ic in range(len(time1)):
1593	tmf1 = tmf1 + time1[ic]
1594	tmf2 = tmf2 + time2[ic]
1595	dtdate1 = dtime.datetime.strptime(tmf1,"%Y-%m-%d_%H:%M:%S")
1596	dtdate2 = dtime.datetime.strptime(tmf2,"%Y-%m-%d_%H:%M:%S")
1597	diffdate12 = dtdate2 - dtdate1
1598	dtime = diffdate12.total_seconds()
1599	print 'dtime:',dtime
1600	else:
1601	print warnmsg
1602	print ' ' + main + ": only 1 time-step for '" + diag + "' !!"
1603	print ' leaving a zero value!'
1604	diagout = var0[:]*0.
1605	dtime=1.
1606
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
1613	ncvar.insert_variable(ncobj, 'pr', diagout/dtime, diagoutd, diagoutvd, newnc)
1614
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
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
1667	# rhs (psfc, t, q) from TimeSeries files
1668	elif diagn == 'TSrhs':
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
1678	diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
1679
1680	ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
1681
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
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
1722	# td (psfc, t, q) from TimeSeries files
1723	elif diagn == 'TStd' or diagn == 'td':
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
1732	diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
1733
1734	ncvar.insert_variable(ncobj, 'tdas', diagout, diagoutd, diagoutvd, newnc)
1735
1736	# td (psfc, t, q) from TimeSeries files
1737	elif diagn == 'TStdC' or diagn == 'tdC':
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
1747	diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
1748
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
1755	ncvar.insert_variable(ncobj, 'tdas', diagout, diagoutd, diagoutvd, newnc)
1756
1757	# wds (u, v)
1758	elif diagn == 'TSwds' or diagn == 'wds' :
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
1766	diagout, diagoutd, diagoutvd = diag.compute_wds(var0,var1,dnamesvar,dvnamesvar)
1767
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
1774	ncvar.insert_variable(ncobj, 'wds', diagout, diagoutd, diagoutvd, newnc)
1775
1776	# wss (u, v)
1777	elif diagn == 'TSwss' or diagn == 'wss':
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
1785	diagout, diagoutd, diagoutvd = diag.compute_wss(var0,var1,dnamesvar,dvnamesvar)
1786
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
1793	ncvar.insert_variable(ncobj, 'wss', diagout, diagoutd, diagoutvd, newnc)
1794
1795	# turbulence (var)
1796	elif diagn == 'turbulence':
1797
1798	var0 = ncobj.variables[depvars][:]
1799
1800	dnamesvar = list(ncobj.variables[depvars].dimensions)
1801	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1802
1803	diagout, diagoutd, diagoutvd = diag.compute_turbulence(var0,dnamesvar,dvnamesvar)
1804	valsvar = gen.variables_values(depvars)
1805
1806	newvarn = depvars + 'turb'
1807	ncvar.insert_variable(ncobj, newvarn, diagout, diagoutd,
1808	diagoutvd, newnc)
1809	varobj = newnc.variables[newvarn]
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
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
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]][:]
1835
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
1845	# WRFbils fom WRF as HFX + LH
1846	elif diagn == 'WRFbils':
1847
1848	var0 = ncobj.variables[depvars[0]][:]
1849	var1 = ncobj.variables[depvars[1]][:]
1850
1851	diagout = var0 + var1
1852	dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1853	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1854
1855	ncvar.insert_variable(ncobj, 'bils', diagout, dnamesvar, dvnamesvar, newnc)
1856
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
1865	var3 = 0.5*(WRFgeop[:,0:dz-1,:,:]+WRFgeop[:,1:dz,:,:])/9.8
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
1888	# WRFclivi WRF water vapour path WRFdens, QICE, QGRAUPEL, QHAIL
1889	elif diagn == 'WRFclivi':
1890
1891	var0 = WRFdens
1892	qtot = ncobj.variables[depvars[1]]
1893	qtotv = qtot[:]
1894	Nspecies = len(depvars) - 2
1895	for iv in range(Nspecies):
1896	if ncobj.variables.has_key(depvars[iv+2]):
1897	var1 = ncobj.variables[depvars[iv+2]][:]
1898	qtotv = qtotv + var1
1899
1900	dnamesvar = list(qtot.dimensions)
1901	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1902
1903	diagout, diagoutd, diagoutvd = diag.compute_clivi(var0, qtotv, dnamesvar,dvnamesvar)
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
1913	# WRFclwvi WRF water cloud-condensed path WRFdens, QCLOUD, QICE, QGRAUPEL, QHAIL
1914	elif diagn == 'WRFclwvi':
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):
1921	if ncobj.variables.has_key(depvars[iv+2]):
1922	var1 = ncobj.variables[depvars[iv+2]]
1923	qtotv = qtotv + var1[:]
1924
1925	dnamesvar = list(qtot.dimensions)
1926	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1927
1928	diagout, diagoutd, diagoutvd = diag.compute_clwvl(var0, qtotv, dnamesvar,dvnamesvar)
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)
1936	ncvar.insert_variable(ncobj, 'clwvi', diagout, diagoutd, diagoutvd, newnc)
1937
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
1964	# WRFgeop geopotential from WRF as PH + PHB
1965	elif diagn == 'WRFgeop':
1966	var0 = ncobj.variables[depvars[0]][:]
1967	var1 = ncobj.variables[depvars[1]][:]
1968
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 = []
1981	diagoutvd = list(dvnames)
1982	for nonvd in NONchkvardims:
1983	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1984	varsadd.append(nonvd)
1985
1986	ncvar.insert_variable(ncobj, 'zg', diagout, dnames, diagoutvd, newnc)
1987
1988	# WRFpotevap_orPM potential evapotranspiration following Penman-Monteith formulation
1989	# implemented in ORCHIDEE (in src_sechiba/enerbil.f90) as: WRFdens, UST, U10, V10, T2, PSFC, QVAPOR
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
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
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
2066	# WRFp pressure from WRF as P + PB
2067	elif diagn == 'WRFp':
2068	var0 = ncobj.variables[depvars[0]][:]
2069	var1 = ncobj.variables[depvars[1]][:]
2070
2071	diagout = var0 + var1
2072	diagoutd = list(ncobj.variables[depvars[0]].dimensions)
2073	diagoutvd = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2074
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)
2080
2081	ncvar.insert_variable(ncobj, 'pres', diagout, diagoutd, diagoutvd, newnc)
2082
2083	# WRFpos
2084	elif diagn == 'WRFpos':
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
2092	elif diagn == 'WRFprw':
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
2100	diagout, diagoutd, diagoutvd = diag.compute_prw(var0, var1, dnamesvar,dvnamesvar)
2101
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)
2108	ncvar.insert_variable(ncobj, 'prw', diagout, diagoutd, diagoutvd, newnc)
2109
2110	# WRFrh (P, T, QVAPOR)
2111	elif diagn == 'WRFrh':
2112
2113	dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
2114	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2115
2116	ncvar.insert_variable(ncobj, 'hur', WRFrh, dnames, dvnames, newnc)
2117
2118	# WRFrhs (PSFC, T2, Q2)
2119	elif diagn == 'WRFrhs':
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
2128	diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
2129	ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
2130
2131	# rvors (u10, v10, WRFpos)
2132	elif diagn == 'WRFrvors':
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
2144	# WRFt (T, P, PB)
2145	elif diagn == 'WRFt':
2146	var0 = ncobj.variables[depvars[0]][:]
2147	var1 = ncobj.variables[depvars[1]][:]
2148	var2 = ncobj.variables[depvars[2]][:]
2149
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
2158	# Removing the nonChecking variable-dimensions from the initial list
2159	varsadd = []
2160	diagoutvd = list(dvnames)
2161	for nonvd in NONchkvardims:
2162	if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2163	varsadd.append(nonvd)
2164
2165	ncvar.insert_variable(ncobj, 'ta', WRFt, dnames, diagoutvd, newnc)
2166
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)
2173	td = ARM3*gammatarh/(ARM2-gammatarh)
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
2187	# WRFtdas (PSFC, T2, Q2)
2188	elif diagn == 'WRFtdas':
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)
2202	tdas = ARM3*gammatarhs/(ARM2-gammatarhs) + 273.15
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
2211	ncvar.insert_variable(ncobj, 'tdas', tdas, dnames, diagoutvd, newnc)
2212
2213	# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
2214	elif diagn == 'WRFua':
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
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
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
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],:])
2234
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
2248	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2249
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)
2256
2257	ncvar.insert_variable(ncobj, 'ua', ua, dnames, diagoutvd, newnc)
2258
2259	# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
2260	elif diagn == 'WRFva':
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
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
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
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
2294	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2295
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)
2303
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
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
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
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],:])
2327
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
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
2356	# WRFtime
2357	elif diagn == 'WRFtime':
2358
2359	diagout = WRFtime
2360
2361	dnamesvar = ['Time']
2362	dvnamesvar = ['Times']
2363
2364	ncvar.insert_variable(ncobj, 'time', diagout, dnamesvar, dvnamesvar, newnc)
2365
2366	# ws (U, V)
2367	elif diagn == 'ws':
2368
2369	# un-staggering variables
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)
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
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
2400	diagout = np.sqrt(unstgvar0unstgvar0 + unstgvar1unstgvar1)
2401
2402	# dnamesvar = ncobj.variables[depvars[0]].dimensions
2403	dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2404
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)
2412
2413	# wss (u10, v10)
2414	elif diagn == 'wss':
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
2426	# WRFheight height from WRF geopotential as WRFGeop/g
2427	elif diagn == 'WRFheight':
2428
2429	diagout = WRFgeop/grav
2430
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)
2437
2438	ncvar.insert_variable(ncobj, 'zhgt', diagout, dnames, diagoutvd, newnc)
2439
2440	# WRFheightrel relative-height from WRF geopotential as WRFgeop(PH + PHB)/g-HGT 'WRFheightrel\|PH@PHB@HGT
2441	elif diagn == 'WRFheightrel':
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):
2449	diagout[:,iz,:,:] = (var0[:,iz,:,:]+ var1[:,iz,:,:])/grav - var2
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
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
2479	# WRFzwind wind extrapolation at a given height using power law computation from WRF
2480	# U, V, WRFz, U10, V10, SINALPHA, COSALPHA, z=[zval]
2481	elif diagn == 'WRFzwind':
2482	var0 = ncobj.variables[depvars[0]][:]
2483	var1 = ncobj.variables[depvars[1]][:]
2484	var2 = WRFz
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,:,:]
2489	var7 = np.float(depvars[7].split('=')[1])
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, \
2500	unvar1, var2, var3, var4, var5, var6, var7, dnames, dvnames)
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
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
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]
2545	# NOTE: only useful for [zval] < 80. m
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
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]
2592	var2 = ncobj.variables[depvars[2]][:]/9.8
2593	var3 = ncobj.getncattr('Station_elevation')
2594
2595	diagout, diagoutd, diagoutvd = diag.Forcompute_zmla_gen(var0,var1,var2,var3, \
2596	dnames,dvnames)
2597
2598	# No need to remove topography height
2599	diagout = diagout + var3
2600
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
2609	else:
2610	print errormsg
2611	print ' ' + main + ": diagnostic '" + diagn + "' not ready!!!"
2612	print ' available diagnostics: ', availdiags
2613	quit(-1)
2614
2615	newnc.sync()
2616	# Adding that additional variables required to compute some diagnostics which
2617	# where not in the original file
2618	print ' adding additional variables...'
2619	for vadd in varsadd:
2620	if not gen.searchInlist(newnc.variables.keys(),vadd) and \
2621	dictcompvars.has_key(vadd):
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)
2635
2636	# end of diagnostics
2637
2638	# Global attributes
2639	##
2640	ncvar.add_global_PyNCplot(newnc, main, None, '2.0')
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 @ 2828

Download in other formats: