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

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

Adding:

  • `1st order 2D horizontal gradient'
File size: 106.0 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
55from optparse import OptionParser
56import numpy as np
57import numpy.ma as ma
58from netCDF4 import Dataset as NetCDFFile
59import os
60import re
61import nc_var_tools as ncvar
62import generic_tools as gen
63import datetime as dtime
64import module_ForDiag as fdin
65import module_ForDef as fdef
66import diag_tools as diag
67
68main = 'diagnostics.py'
69errormsg = 'ERROR -- error -- ERROR -- error'
70warnmsg = 'WARNING -- warning -- WARNING -- warning'
71
72# Constants
73grav = 9.81
74
75
76####### ###### ##### #### ### ## #
77comboinf="\nIF -d 'variable_combo', provides information of the combination to obtain -v [varn] with the ASCII file with the combinations as -f [combofile]"
78
79parser = OptionParser()
80parser.add_option("-f", "--netCDF_file", dest="ncfile", help="file to use", metavar="FILE")
81parser.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")
84parser.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    #######
92availdiags = ['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
102methods = ['accum', 'deaccum']
103
104# Variables not to check
105NONcheckingvars = ['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
117NONcheckdepvars = ['accum', 'deaccum', 'timeoverthres', 'WRF_denszint',              \
118  'WRFzwind_log', 'WRFzwind', 'WRFzwindMO']
119
120NONchkvardims = ['WRFtime']
121
122ofile = 'diagnostics.nc'
123
124dimns = opts.dimns
125varns = opts.varns
126
127# Special method. knowing variable combination
128##
129if 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
136if 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
142if 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
148if 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
154if not os.path.isfile(opts.ncfile):
155    print errormsg
156    print '   ' + main + ": file '" + opts.ncfile + "' does not exist !!"
157    quit(-1)
158
159ncobj = NetCDFFile(opts.ncfile, 'r')
160
161# Looking for specific variables that might be use in more than one diagnostic
162WRFgeop_compute = False
163WRFp_compute = False
164WRFt_compute = False
165WRFrh_compute = False
166WRFght_compute = False
167WRFdens_compute = False
168WRFpos_compute = False
169WRFtime_compute = False
170WRFz_compute = False
171WRFdx_compute = False
172WRFdy_compute = False
173WRFdz_compute = False
174WRFdxdy_compute = False
175WRFdxdywps_compute = False
176LONLATdxdy_compute = False
177LONLATdx_compute = False
178LONLATdy_compute = False
179UNua_compute = False
180UNva_compute = False
181UNwa_compute = False
182
183# File creation
184newnc = NetCDFFile(ofile,'w')
185
186# dimensions
187dimvalues = dimns.split(',')
188dnames = []
189dvnames = []
190
191for 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
225diags = varns.split(',')
226Ndiags = len(diags)
227
228for 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
293dictcompvars = {}
294if 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
304if 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
314if 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
323if 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.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
332    data2 = 0.622*data1/(0.01*p-(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
341if 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
353if 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 = dxval*dyval*mapfac
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
381if WRFpos_compute:
382# WRF positions from the lowest-leftest corner of the matrix
383    print '    ' + main + ': computing position from MAPFAC_M as sqrt(DY*j**2 + ' +  \
384      'DX*x**2)*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((disty*j)**2. + (distx*i)**2.)/mapfac[0,j,i]
406#            WRFpos[0,j,i] = np.sqrt((disty*j)**2. + (distx*i)**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
412if 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
453if 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
472if 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(dx**2+dy**2)'
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(WRFdx**2 + WRFdy**2)
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
504if 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(dx**2+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(WRFdx**2 + WRFdy**2)
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
536if 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
559if 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])*raddeg*Rearth*' +    \
562      'cos(abs(lat[i])); dy=(lat[j+1]-lat[i])*raddeg*Rearth; ds=sqrt(dx**2+dy**2); '+\
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])*raddeg*Rearth*           \
581      np.cos(np.abs(WRFlat[:,0:dimx-1]*raddeg))
582    WRFdy[0:dimy-1,:]=(WRFlat[1:dimy,:]-WRFlat[0:dimy-1,:])*raddeg*Rearth
583    WRFds = np.sqrt(WRFdx**2 + WRFdy**2)
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
596if 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
617if 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
638if 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### ## #
662print '  ' + main + ' ...'
663varsadd = []
664
665Varns = ncobj.variables.keys()
666Varns.sort()
667
668for 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 = ncvar.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        var1 = ncobj.variables[depvars[1]]
1189        var2 = ncobj.variables[depvars[2]]
1190
1191        dnamesvar = var0.dimensions
1192        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1193
1194        diagout, diagoutd, diagoutvd = diag.Forcompute_gradient2Dh(var0, var1, var2, \
1195          dnamesvar, dvnamesvar)
1196        # Removing the nonChecking variable-dimensions from the initial list
1197        varsadd = []
1198        diagoutvd = list(dvnames)
1199        for nonvd in NONchkvardims:
1200            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1201            varsadd.append(nonvd)
1202
1203        CFvarn = ncvar.variables_values(depvars[0])[0]
1204        ncvar.insert_variable(ncobj, depvars[0], diagout, diagoutd, diagoutvd,       \
1205          newnc, gen.fillValueF)
1206        newnc.sync()
1207        ovar = newnc.variables[depvars[0]]
1208        varu = ovar.units
1209        newnc.renameVariable(depvars[0], depvars[0]+'2Dhgrad')
1210        ovar.setncattr('units', varu+'ds-1')
1211
1212# LMDZrh (pres, t, r)
1213    elif diagn == 'LMDZrh':
1214           
1215        var0 = ncobj.variables[depvars[0]][:]
1216        var1 = ncobj.variables[depvars[1]][:]
1217        var2 = ncobj.variables[depvars[2]][:]
1218
1219        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnames,dvnames)
1220        ncvar.insert_variable(ncobj, 'hur', diagout, diagoutd, diagoutvd, newnc)
1221
1222# LMDZrhs (psol, t2m, q2m)
1223    elif diagn == 'LMDZrhs':
1224           
1225        var0 = ncobj.variables[depvars[0]][:]
1226        var1 = ncobj.variables[depvars[1]][:]
1227        var2 = ncobj.variables[depvars[2]][:]
1228
1229        dnamesvar = ncobj.variables[depvars[0]].dimensions
1230        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1231
1232        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
1233
1234        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
1235
1236# range_faces: LON, LAT, HGT, WRFdxdy, 'face:['WE'/'SN']:[dsfilt]:[dsnewrange]:[hvalleyrange]'
1237    elif diagn == 'range_faces':
1238           
1239        var0 = ncobj.variables[depvars[0]][:]
1240        var1 = ncobj.variables[depvars[1]][:]
1241        var2 = ncobj.variables[depvars[2]][:]
1242        face = depvars[4].split(':')[1]
1243        dsfilt = np.float(depvars[4].split(':')[2])
1244        dsnewrange = np.float(depvars[4].split(':')[3])
1245        hvalleyrange = np.float(depvars[4].split(':')[4])
1246
1247        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
1248        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1249        lon, lat = gen.lonlat2D(var0, var1)
1250        if len(var2.shape) == 3:
1251            print warnmsg
1252            print '  ' + diagn + ": shapping to 2D variable '" + depvars[2] + "' !!"
1253            hgt = var2[0,:,:]
1254            dnamesvar.pop(0)
1255            dvnamesvar.pop(0)           
1256        else:
1257            hgt = var2[:]
1258
1259        orogmax, ptorogmax, dhgt, peaks, valleys, origfaces, diagout, diagoutd,      \
1260          diagoutvd, rng, rngorogmax, ptrngorogmax= diag.Forcompute_range_faces(lon, \
1261           lat, hgt, WRFdx, WRFdy, WRFds, face, dsfilt, dsnewrange, hvalleyrange,    \
1262           dnamesvar, dvnamesvar)
1263
1264        # Removing the nonChecking variable-dimensions from the initial list
1265        varsadd = []
1266        diagoutvd = list(dvnames)
1267        for nonvd in NONchkvardims:
1268            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1269            varsadd.append(nonvd)
1270
1271        ncvar.insert_variable(ncobj, 'dx', WRFdx, diagoutd, diagoutvd, newnc)
1272        ncvar.insert_variable(ncobj, 'dy', WRFdy, diagoutd, diagoutvd, newnc)
1273        ncvar.insert_variable(ncobj, 'ds', WRFds, diagoutd, diagoutvd, newnc)
1274
1275        # adding variables to output file
1276        if face == 'WE': axis = 'lon'
1277        elif face == 'SN': axis = 'lat'
1278
1279        ncvar.insert_variable(ncobj, 'range', rng, diagoutd, diagoutvd, newnc,       \
1280          fill=gen.fillValueI)
1281        ovar = newnc.variables['range']
1282        ncvar.set_attribute(ovar, 'deriv', axis)
1283
1284        ncvar.insert_variable(ncobj, 'orogmax', rngorogmax, diagoutd, diagoutvd,     \
1285          newnc, fill=gen.fillValueF)
1286        newnc.renameVariable('orogmax', 'rangeorogmax')
1287        ovar = newnc.variables['rangeorogmax']
1288        ncvar.set_attribute(ovar, 'deriv', axis)
1289        stdn = ovar.standard_name
1290        ncvar.set_attribute(ovar, 'standard_name', 'range_' + stdn)
1291        Ln = ovar.long_name
1292        ncvar.set_attribute(ovar, 'long_name', 'range ' + stdn)
1293
1294        ncvar.insert_variable(ncobj, 'ptorogmax', ptrngorogmax, diagoutd, diagoutvd, \
1295          newnc)
1296        newnc.renameVariable('ptorogmax', 'rangeptorogmax')
1297        ovar = newnc.variables['rangeptorogmax']
1298        ncvar.set_attribute(ovar, 'deriv', axis)
1299        stdn = ovar.standard_name
1300        ncvar.set_attribute(ovar, 'standard_name', 'range_' + stdn)
1301        Ln = ovar.long_name
1302        ncvar.set_attribute(ovar, 'long_name', 'range ' + stdn)
1303
1304        ncvar.insert_variable(ncobj, 'orogmax', orogmax, diagoutd, diagoutvd,        \
1305          newnc)
1306        ovar = newnc.variables['orogmax']
1307        ncvar.set_attribute(ovar, 'deriv', axis)
1308
1309        ncvar.insert_variable(ncobj, 'ptorogmax', ptorogmax, diagoutd, diagoutvd,    \
1310          newnc)
1311        ovar = newnc.variables['ptorogmax']
1312        ncvar.set_attribute(ovar, 'deriv', axis)
1313
1314        ncvar.insert_variable(ncobj, 'orogderiv', dhgt, diagoutd, diagoutvd, newnc)
1315        ovar = newnc.variables['orogderiv']
1316        ncvar.set_attribute(ovar, 'deriv', axis)
1317
1318        ncvar.insert_variable(ncobj, 'peak', peaks, diagoutd, diagoutvd, newnc)
1319        ncvar.insert_variable(ncobj, 'valley', valleys, diagoutd, diagoutvd, newnc)
1320
1321        ncvar.insert_variable(ncobj, 'rangefaces', diagout, diagoutd, diagoutvd,    \
1322          newnc)
1323        newnc.renameVariable('rangefaces', 'rangefacesfilt')
1324        ovar = newnc.variables['rangefacesfilt']
1325        ncvar.set_attribute(ovar, 'face', face)
1326        ncvar.set_attributek(ovar, 'dist_filter', dsfilt, 'F')
1327
1328        ncvar.insert_variable(ncobj, 'rangefaces', origfaces, diagoutd, diagoutvd,  \
1329          newnc, fill=gen.fillValueI)
1330        ovar = newnc.variables['rangefaces']
1331        ncvar.set_attribute(ovar, 'face', face)
1332        ncvar.set_attributek(ovar, 'dist_newrange', dsnewrange, 'F')
1333        ncvar.set_attributek(ovar, 'h_valley_newrange', hvalleyrange, 'F')
1334
1335# cell_bnds: grid cell bounds from lon, lat from a reglar lon/lat projection  as
1336#   intersection of their related parallels and meridians
1337    elif diagn == 'reglonlatbnds':
1338           
1339        var00 = ncobj.variables[depvars[0]][:]
1340        var01 = ncobj.variables[depvars[1]][:]
1341
1342        var0, var1 = gen.lonlat2D(var00,var01)
1343
1344        dnamesvar = []
1345        dnamesvar.append('bnds')
1346        if (len(var00.shape) == 3):
1347            dnamesvar.append(ncobj.variables[depvars[0]].dimensions[1])
1348            dnamesvar.append(ncobj.variables[depvars[0]].dimensions[2])
1349        elif (len(var00.shape) == 2):
1350            dnamesvar.append(ncobj.variables[depvars[0]].dimensions[0])
1351            dnamesvar.append(ncobj.variables[depvars[0]].dimensions[1])
1352        elif (len(var00.shape) == 1):
1353            dnamesvar.append(ncobj.variables[depvars[0]].dimensions[0])
1354            dnamesvar.append(ncobj.variables[depvars[1]].dimensions[0])
1355        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1356
1357        cellbndsx, cellbndsy, diagoutd, diagoutvd = diag.Forcompute_cellbndsreg(var0,\
1358          var1, dnamesvar, dvnamesvar)
1359
1360        # Removing the nonChecking variable-dimensions from the initial list
1361        varsadd = []
1362        diagoutvd = list(dvnames)
1363        for nonvd in NONchkvardims:
1364            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1365            varsadd.append(nonvd)
1366        # creation of bounds dimension
1367        newdim = newnc.createDimension('bnds', 4)
1368
1369        ncvar.insert_variable(ncobj, 'lon_bnds', cellbndsx, diagoutd, diagoutvd, newnc)
1370        ncvar.insert_variable(ncobj, 'lat_bnds', cellbndsy, diagoutd, diagoutvd, newnc)
1371
1372# tws: Bet Wulb temperature following Stull 2011 (tas, hurs)
1373    elif diagn == 'tws':
1374           
1375        var0 = ncobj.variables[depvars[0]][:]
1376        var1 = ncobj.variables[depvars[1]][:]
1377
1378        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1379        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1380        diagoutd = dnames
1381        diagoutvd = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1382
1383        diagout = diag.var_tws_S11(var0,var1)
1384        ncvar.insert_variable(ncobj, 'tws', diagout, diagoutd, diagoutvd, newnc)
1385
1386# cell_bnds: grid cell bounds from XLONG_U, XLAT_U, XLONG_V, XLAT_V as intersection
1387#   of their related parallels and meridians
1388    elif diagn == 'WRFbnds':
1389           
1390        var0 = ncobj.variables[depvars[0]][0,:,:]
1391        var1 = ncobj.variables[depvars[1]][0,:,:]
1392        var2 = ncobj.variables[depvars[2]][0,:,:]
1393        var3 = ncobj.variables[depvars[3]][0,:,:]
1394
1395        dnamesvar = []
1396        dnamesvar.append('bnds')
1397        dnamesvar.append(ncobj.variables[depvars[0]].dimensions[1])
1398        dnamesvar.append(ncobj.variables[depvars[2]].dimensions[2])
1399        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1400
1401        cellbndsx, cellbndsy, diagoutd, diagoutvd = diag.Forcompute_cellbnds(var0,   \
1402          var1, var2, var3, dnamesvar, dvnamesvar)
1403
1404        # Removing the nonChecking variable-dimensions from the initial list
1405        varsadd = []
1406        diagoutvd = list(dvnames)
1407        for nonvd in NONchkvardims:
1408            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1409            varsadd.append(nonvd)
1410        # creation of bounds dimension
1411        newdim = newnc.createDimension('bnds', 4)
1412
1413        ncvar.insert_variable(ncobj, 'lon_bnds', cellbndsx, diagoutd, diagoutvd, newnc)
1414        newnc.sync()
1415        ncvar.insert_variable(ncobj, 'lat_bnds', cellbndsy, diagoutd, diagoutvd, newnc)
1416        newnc.sync()
1417
1418        if newnc.variables.has_key('XLONG'):
1419            ovar = newnc.variables['XLONG']
1420            ovar.setncattr('bounds', 'lon_bnds lat_bnds')
1421            ovar = newnc.variables['XLAT']
1422            ovar.setncattr('bounds', 'lon_bnds lat_bnds')
1423        elif newnc.variables.has_key('XLONG_M'):
1424            ovar = newnc.variables['XLONG_M']
1425            ovar.setncattr('bounds', 'lon_bnds lat_bnds')
1426            ovar = newnc.variables['XLAT_M']
1427            ovar.setncattr('bounds', 'lon_bnds lat_bnds')
1428        else:
1429            print errormsg
1430            print '  ' + fname + ": error computing diagnostic '" + diagn + "' !!"
1431            print "    neither: 'XLONG', 'XLONG_M' have been found"
1432            quit(-1)
1433
1434# mrso: total soil moisture SMOIS, DZS
1435    elif diagn == 'WRFmrso':
1436           
1437        var0 = ncobj.variables[depvars[0]][:]
1438        var10 = ncobj.variables[depvars[1]][:]
1439        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1440        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1441
1442        var1 = var0.copy()*0.
1443        var2 = var0.copy()*0.+1.
1444        # Must be a better way....
1445        for j in range(var0.shape[2]):
1446          for i in range(var0.shape[3]):
1447              var1[:,:,j,i] = var10
1448
1449        diagout, diagoutd, diagoutvd = diag.Forcompute_zint(var0, var1, var2,        \
1450          dnamesvar, dvnamesvar)
1451
1452        # Removing the nonChecking variable-dimensions from the initial list
1453        varsadd = []
1454        diagoutvd = list(dvnames)
1455        for nonvd in NONchkvardims:
1456            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1457            varsadd.append(nonvd)
1458        ncvar.insert_variable(ncobj, 'mrso', diagout, diagoutd, diagoutvd, newnc)
1459
1460# mrsos: First layer soil moisture SMOIS, DZS
1461    elif diagn == 'WRFmrsos':
1462           
1463        var0 = ncobj.variables[depvars[0]][:]
1464        var1 = ncobj.variables[depvars[1]][:]
1465        diagoutd = list(ncobj.variables[depvars[0]].dimensions)
1466        diagoutvd = ncvar.var_dim_dimv(diagoutd,dnames,dvnames)
1467
1468        diagoutd.pop(1)
1469        diagoutvd.pop(1)
1470
1471        diagout= np.zeros((var0.shape[0],var0.shape[2],var0.shape[3]), dtype=np.float)
1472
1473        # Must be a better way....
1474        for j in range(var0.shape[2]):
1475          for i in range(var0.shape[3]):
1476              diagout[:,j,i] = var0[:,0,j,i]*var1[:,0]
1477
1478        # Removing the nonChecking variable-dimensions from the initial list
1479        varsadd = []
1480        diagoutvd = list(dvnames)
1481        for nonvd in NONchkvardims:
1482            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1483            varsadd.append(nonvd)
1484        ncvar.insert_variable(ncobj, 'mrsos', diagout, diagoutd, diagoutvd, newnc)
1485
1486# mslp: mean sea level pressure (pres, psfc, terrain, temp, qv)
1487    elif diagn == 'mslp' or diagn == 'WRFmslp':
1488           
1489        var1 = ncobj.variables[depvars[1]][:]
1490        var2 = ncobj.variables[depvars[2]][:]
1491        var4 = ncobj.variables[depvars[4]][:]
1492
1493        if diagn == 'WRFmslp':
1494            var0 = WRFp
1495            var3 = WRFt
1496            dnamesvar = ncobj.variables['P'].dimensions
1497        else:
1498            var0 = ncobj.variables[depvars[0]][:]
1499            var3 = ncobj.variables[depvars[3]][:]
1500            dnamesvar = ncobj.variables[depvars[0]].dimensions
1501
1502        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1503
1504        diagout, diagoutd, diagoutvd = diag.compute_mslp(var0, var1, var2, var3, var4,    \
1505          dnamesvar, dvnamesvar)
1506
1507        # Removing the nonChecking variable-dimensions from the initial list
1508        varsadd = []
1509        diagoutvd = list(dvnames)
1510        for nonvd in NONchkvardims:
1511            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1512            varsadd.append(nonvd)
1513        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
1514
1515# WRFtws: Bet Wulb temperature following Stull 2011 (PSFC, T2, Q2)
1516    elif diagn == 'WRFtws':
1517
1518        var0 = ncobj.variables[depvars[0]][:]
1519        var1 = ncobj.variables[depvars[1]][:]
1520        var2 = ncobj.variables[depvars[2]][:]
1521
1522        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
1523        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1524
1525        hurs, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
1526           
1527        diagout = diag.var_tws_S11(var1, hurs)
1528        # Removing the nonChecking variable-dimensions from the initial list
1529        varsadd = []
1530        diagoutvd = list(dvnames)
1531        for nonvd in NONchkvardims:
1532            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1533            varsadd.append(nonvd)
1534
1535        ncvar.insert_variable(ncobj, 'tws', diagout, diagoutd, diagoutvd, newnc)
1536
1537# OMEGAw (omega, p, t) from NCL formulation (https://www.ncl.ucar.edu/Document/Functions/Contributed/omega_to_w.shtml)
1538    elif diagn == 'OMEGAw':
1539           
1540        var0 = ncobj.variables[depvars[0]][:]
1541        var1 = ncobj.variables[depvars[1]][:]
1542        var2 = ncobj.variables[depvars[2]][:]
1543
1544        dnamesvar = ncobj.variables[depvars[0]].dimensions
1545        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1546
1547        diagout, diagoutd, diagoutvd = diag.compute_OMEGAw(var0,var1,var2,dnamesvar,dvnamesvar)
1548
1549        ncvar.insert_variable(ncobj, 'wa', diagout, diagoutd, diagoutvd, newnc)
1550
1551# raintot: instantaneous total precipitation from WRF as (RAINC + RAINC + RAINSH) / dTime
1552    elif diagn == 'RAINTOT':
1553
1554        var0 = ncobj.variables[depvars[0]]
1555        var1 = ncobj.variables[depvars[1]]
1556        var2 = ncobj.variables[depvars[2]]
1557
1558        if depvars[3] != 'WRFtime':
1559            var3 = ncobj.variables[depvars[3]]
1560        else:
1561            var3 = np.arange(var0.shape[0], dtype=int)
1562
1563        var = var0[:] + var1[:] + var2[:]
1564
1565        dnamesvar = var0.dimensions
1566        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1567
1568        diagout, diagoutd, diagoutvd = diag.compute_deaccum(var,dnamesvar,dvnamesvar)
1569
1570# Transforming to a flux
1571        if var3.shape[0] > 1:
1572            if depvars[3] != 'WRFtime':
1573                dtimeunits = var3.getncattr('units')
1574                tunits = dtimeunits.split(' ')[0]
1575   
1576                dtime = (var3[1] - var3[0])*diag.timeunits_seconds(tunits)
1577            else:
1578                var3 = ncobj.variables['Times']
1579                time1 = var3[0,:]
1580                time2 = var3[1,:]
1581                tmf1 = ''
1582                tmf2 = ''
1583                for ic in range(len(time1)):
1584                    tmf1 = tmf1 + time1[ic]
1585                    tmf2 = tmf2 + time2[ic]
1586                dtdate1 = dtime.datetime.strptime(tmf1,"%Y-%m-%d_%H:%M:%S")
1587                dtdate2 = dtime.datetime.strptime(tmf2,"%Y-%m-%d_%H:%M:%S")
1588                diffdate12 = dtdate2 - dtdate1
1589                dtime = diffdate12.total_seconds()
1590                print 'dtime:',dtime
1591        else:
1592            print warnmsg
1593            print '  ' + main + ": only 1 time-step for '" + diag + "' !!"
1594            print '    leaving a zero value!'
1595            diagout = var0[:]*0.
1596            dtime=1.
1597
1598        # Removing the nonChecking variable-dimensions from the initial list
1599        varsadd = []
1600        for nonvd in NONchkvardims:
1601            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
1602            varsadd.append(nonvd)
1603           
1604        ncvar.insert_variable(ncobj, 'pr', diagout/dtime, diagoutd, diagoutvd, newnc)
1605
1606# timemax ([varname], time). When a given variable [varname] got its maximum
1607    elif diagn == 'timemax':
1608           
1609        var0 = ncobj.variables[depvars[0]][:]
1610        var1 = ncobj.variables[depvars[1]][:]
1611
1612        otime = ncobj.variables[depvars[1]]
1613
1614        dnamesvar = ncobj.variables[depvars[0]].dimensions
1615        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1616
1617        diagout, diagoutd, diagoutvd  = diag.var_timemax(var0, var1, dnames,         \
1618          dvnames)
1619
1620        ncvar.insert_variable(ncobj, 'timemax', diagout, diagoutd, diagoutvd, newnc, \
1621          fill=gen.fillValueF)
1622        # Getting the right units
1623        ovar = newnc.variables['timemax']
1624        if gen.searchInlist(otime.ncattrs(), 'units'): 
1625            tunits = otime.getncattr('units')
1626            ncvar.set_attribute(ovar, 'units', tunits)
1627            newnc.sync()
1628        ncvar.set_attribute(ovar, 'variable', depvars[0])
1629
1630# timeoverthres ([varname], time, [value], [CFvarn]). When a given variable [varname]   
1631#   overpass a given [value]. Being [CFvarn] the name of the diagnostics in
1632#   variables_values.dat
1633    elif diagn == 'timeoverthres':
1634           
1635        var0 = ncobj.variables[depvars[0]][:]
1636        var1 = ncobj.variables[depvars[1]][:]
1637        var2 = np.float(depvars[2])
1638        var3 = depvars[3]
1639
1640        otime = ncobj.variables[depvars[1]]
1641
1642        dnamesvar = ncobj.variables[depvars[0]].dimensions
1643        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1644
1645        diagout, diagoutd, diagoutvd  = diag.var_timeoverthres(var0, var1, var2,     \
1646          dnames, dvnames)
1647
1648        ncvar.insert_variable(ncobj, var3, diagout, diagoutd, diagoutvd, newnc, \
1649          fill=gen.fillValueF)
1650        # Getting the right units
1651        ovar = newnc.variables[var3]
1652        if gen.searchInlist(otime.ncattrs(), 'units'): 
1653            tunits = otime.getncattr('units')
1654            ncvar.set_attribute(ovar, 'units', tunits)
1655            newnc.sync()
1656        ncvar.set_attribute(ovar, 'overpassed_threshold', var2)
1657
1658# rhs (psfc, t, q) from TimeSeries files
1659    elif diagn == 'TSrhs':
1660           
1661        p0=100000.
1662        var0 = ncobj.variables[depvars[0]][:]
1663        var1 = (ncobj.variables[depvars[1]][:])*(var0/p0)**(2./7.)
1664        var2 = ncobj.variables[depvars[2]][:]
1665
1666        dnamesvar = ncobj.variables[depvars[0]].dimensions
1667        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1668
1669        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
1670
1671        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
1672
1673# rhs (psfc, t, q) from tas, tds
1674    elif diagn == 'rhs_tas_tds':
1675           
1676        var0 = ncobj.variables[depvars[0]][:]
1677        var1 = ncobj.variables[depvars[1]][:]
1678
1679        dnamesvar = ncobj.variables[depvars[0]].dimensions
1680        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1681
1682        diagout, diagoutd, diagoutvd = diag.var_hur_tas_tds(var0,var1,dnamesvar,     \
1683          dvnamesvar)
1684
1685        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
1686
1687# slw: total soil liquid water SH2O, DZS
1688    elif diagn == 'WRFslw':
1689           
1690        var0 = ncobj.variables[depvars[0]][:]
1691        var10 = ncobj.variables[depvars[1]][:]
1692        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1693        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1694
1695        var1 = var0.copy()*0.
1696        var2 = var0.copy()*0.+1.
1697        # Must be a better way....
1698        for j in range(var0.shape[2]):
1699          for i in range(var0.shape[3]):
1700              var1[:,:,j,i] = var10
1701
1702        diagout, diagoutd, diagoutvd = diag.Forcompute_zint(var0, var1, var2,        \
1703          dnamesvar, dvnamesvar)
1704
1705        # Removing the nonChecking variable-dimensions from the initial list
1706        varsadd = []
1707        diagoutvd = list(dvnames)
1708        for nonvd in NONchkvardims:
1709            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1710            varsadd.append(nonvd)
1711        ncvar.insert_variable(ncobj, 'slw', diagout, diagoutd, diagoutvd, newnc)
1712
1713# td (psfc, t, q) from TimeSeries files
1714    elif diagn == 'TStd' or diagn == 'td':
1715           
1716        var0 = ncobj.variables[depvars[0]][:]
1717        var1 = ncobj.variables[depvars[1]][:] - 273.15
1718        var2 = ncobj.variables[depvars[2]][:]
1719
1720        dnamesvar = ncobj.variables[depvars[0]].dimensions
1721        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1722
1723        diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
1724
1725        ncvar.insert_variable(ncobj, 'tdas', diagout, diagoutd, diagoutvd, newnc)
1726
1727# td (psfc, t, q) from TimeSeries files
1728    elif diagn == 'TStdC' or diagn == 'tdC':
1729           
1730        var0 = ncobj.variables[depvars[0]][:]
1731# Temperature is already in degrees Celsius
1732        var1 = ncobj.variables[depvars[1]][:]
1733        var2 = ncobj.variables[depvars[2]][:]
1734
1735        dnamesvar = ncobj.variables[depvars[0]].dimensions
1736        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1737
1738        diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
1739
1740        # Removing the nonChecking variable-dimensions from the initial list
1741        varsadd = []
1742        for nonvd in NONchkvardims:
1743            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1744            varsadd.append(nonvd)
1745
1746        ncvar.insert_variable(ncobj, 'tdas', diagout, diagoutd, diagoutvd, newnc)
1747
1748# wds (u, v)
1749    elif diagn == 'TSwds' or diagn == 'wds' :
1750 
1751        var0 = ncobj.variables[depvars[0]][:]
1752        var1 = ncobj.variables[depvars[1]][:]
1753
1754        dnamesvar = ncobj.variables[depvars[0]].dimensions
1755        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1756
1757        diagout, diagoutd, diagoutvd = diag.compute_wds(var0,var1,dnamesvar,dvnamesvar)
1758
1759        # Removing the nonChecking variable-dimensions from the initial list
1760        varsadd = []
1761        for nonvd in NONchkvardims:
1762            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1763            varsadd.append(nonvd)
1764
1765        ncvar.insert_variable(ncobj, 'wds', diagout, diagoutd, diagoutvd, newnc)
1766
1767# wss (u, v)
1768    elif diagn == 'TSwss' or diagn == 'wss':
1769           
1770        var0 = ncobj.variables[depvars[0]][:]
1771        var1 = ncobj.variables[depvars[1]][:]
1772
1773        dnamesvar = ncobj.variables[depvars[0]].dimensions
1774        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1775
1776        diagout, diagoutd, diagoutvd = diag.compute_wss(var0,var1,dnamesvar,dvnamesvar)
1777
1778        # Removing the nonChecking variable-dimensions from the initial list
1779        varsadd = []
1780        for nonvd in NONchkvardims:
1781            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1782            varsadd.append(nonvd)
1783
1784        ncvar.insert_variable(ncobj, 'wss', diagout, diagoutd, diagoutvd, newnc)
1785
1786# turbulence (var)
1787    elif diagn == 'turbulence':
1788
1789        var0 = ncobj.variables[depvars][:]
1790
1791        dnamesvar = list(ncobj.variables[depvars].dimensions)
1792        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1793
1794        diagout, diagoutd, diagoutvd = diag.compute_turbulence(var0,dnamesvar,dvnamesvar)
1795        valsvar = gen.variables_values(depvars)
1796
1797        newvarn = depvars + 'turb'
1798        ncvar.insert_variable(ncobj, newvarn, diagout, diagoutd, 
1799          diagoutvd, newnc)
1800        varobj = newnc.variables[newvarn]
1801        attrv = varobj.long_name
1802        attr = varobj.delncattr('long_name')
1803        newattr = ncvar.set_attribute(varobj, 'long_name', attrv +                   \
1804          " Taylor decomposition turbulence term")
1805
1806# ua va from ws wd (deg)
1807    elif diagn == 'uavaFROMwswd':
1808           
1809        var0 = ncobj.variables[depvars[0]][:]
1810        var1 = ncobj.variables[depvars[1]][:]
1811
1812        ua = var0*np.cos(var1*np.pi/180.)
1813        va = var0*np.sin(var1*np.pi/180.)
1814
1815        dnamesvar = ncobj.variables[depvars[0]].dimensions
1816        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1817
1818        ncvar.insert_variable(ncobj, 'ua', ua, dnamesvar, dvnamesvar, newnc)
1819        ncvar.insert_variable(ncobj, 'va', va, dnamesvar, dvnamesvar, newnc)
1820
1821# ua va from obs ws wd (deg)
1822    elif diagn == 'uavaFROMobswswd':
1823           
1824        var0 = ncobj.variables[depvars[0]][:]
1825        var1 = ncobj.variables[depvars[1]][:]
1826
1827        ua = var0*np.cos((var1+180.)*np.pi/180.)
1828        va = var0*np.sin((var1+180.)*np.pi/180.)
1829
1830        dnamesvar = ncobj.variables[depvars[0]].dimensions
1831        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1832
1833        ncvar.insert_variable(ncobj, 'ua', ua, dnamesvar, dvnamesvar, newnc)
1834        ncvar.insert_variable(ncobj, 'va', va, dnamesvar, dvnamesvar, newnc)
1835
1836# WRFbils fom WRF as HFX + LH
1837    elif diagn == 'WRFbils':
1838           
1839        var0 = ncobj.variables[depvars[0]][:]
1840        var1 = ncobj.variables[depvars[1]][:]
1841
1842        diagout = var0 + var1
1843        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1844        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1845
1846        ncvar.insert_variable(ncobj, 'bils', diagout, dnamesvar, dvnamesvar, newnc)
1847
1848# WRFcape_afwa CAPE, CIN, ZLFC, PLFC, LI following WRF 'phys/module_diaf_afwa.F'
1849#   methodology as WRFt, WRFrh, WRFp, WRFgeop, HGT
1850    elif diagn == 'WRFcape_afwa':
1851        var0 = WRFt
1852        var1 = WRFrh
1853        var2 = WRFp
1854        dz = WRFgeop.shape[1]
1855        # de-staggering
1856        var3 = 0.5*(WRFgeop[:,0:dz-1,:,:]+WRFgeop[:,1:dz,:,:])/9.8
1857        var4 = ncobj.variables[depvars[4]][0,:,:]
1858
1859        dnamesvar = list(ncobj.variables['T'].dimensions)
1860        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1861
1862        diagout = np.zeros(var0.shape, dtype=np.float)
1863        diagout1, diagout2, diagout3, diagout4, diagout5, diagoutd, diagoutvd =      \
1864          diag.Forcompute_cape_afwa(var0, var1, var2, var3, var4, 3, dnamesvar,      \
1865          dvnamesvar)
1866
1867        # Removing the nonChecking variable-dimensions from the initial list
1868        varsadd = []
1869        for nonvd in NONchkvardims:
1870            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1871            varsadd.append(nonvd)
1872
1873        ncvar.insert_variable(ncobj, 'cape', diagout1, diagoutd, diagoutvd, newnc)
1874        ncvar.insert_variable(ncobj, 'cin', diagout2, diagoutd, diagoutvd, newnc)
1875        ncvar.insert_variable(ncobj, 'zlfc', diagout3, diagoutd, diagoutvd, newnc)
1876        ncvar.insert_variable(ncobj, 'plfc', diagout4, diagoutd, diagoutvd, newnc)
1877        ncvar.insert_variable(ncobj, 'li', diagout5, diagoutd, diagoutvd, newnc)
1878
1879# WRFclivi WRF water vapour path WRFdens, QICE, QGRAUPEL, QHAIL
1880    elif diagn == 'WRFclivi':
1881           
1882        var0 = WRFdens
1883        qtot = ncobj.variables[depvars[1]]
1884        qtotv = qtot[:]
1885        Nspecies = len(depvars) - 2
1886        for iv in range(Nspecies):
1887            if ncobj.variables.has_key(depvars[iv+2]):
1888                var1 = ncobj.variables[depvars[iv+2]][:]
1889                qtotv = qtotv + var1
1890
1891        dnamesvar = list(qtot.dimensions)
1892        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1893
1894        diagout, diagoutd, diagoutvd = diag.compute_clivi(var0, qtotv, dnamesvar,dvnamesvar)
1895
1896        # Removing the nonChecking variable-dimensions from the initial list
1897        varsadd = []
1898        diagoutvd = list(dvnames)
1899        for nonvd in NONchkvardims:
1900            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1901            varsadd.append(nonvd)
1902        ncvar.insert_variable(ncobj, 'clivi', diagout, diagoutd, diagoutvd, newnc)
1903
1904# WRFclwvi WRF water cloud-condensed path WRFdens, QCLOUD, QICE, QGRAUPEL, QHAIL
1905    elif diagn == 'WRFclwvi':
1906           
1907        var0 = WRFdens
1908        qtot = ncobj.variables[depvars[1]]
1909        qtotv = ncobj.variables[depvars[1]]
1910        Nspecies = len(depvars) - 2
1911        for iv in range(Nspecies):
1912            if ncobj.variables.has_key(depvars[iv+2]):
1913                var1 = ncobj.variables[depvars[iv+2]]
1914                qtotv = qtotv + var1[:]
1915
1916        dnamesvar = list(qtot.dimensions)
1917        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1918
1919        diagout, diagoutd, diagoutvd = diag.compute_clwvl(var0, qtotv, dnamesvar,dvnamesvar)
1920
1921        # Removing the nonChecking variable-dimensions from the initial list
1922        varsadd = []
1923        diagoutvd = list(dvnames)
1924        for nonvd in NONchkvardims:
1925            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1926            varsadd.append(nonvd)
1927        ncvar.insert_variable(ncobj, 'clwvi', diagout, diagoutd, diagoutvd, newnc)
1928
1929# WRF_denszint WRF vertical integration as WRFdens, sum(Q[water species1], ..., Q[water speciesN]), varn=[varN]
1930    elif diagn == 'WRF_denszint':
1931           
1932        var0 = WRFdens
1933        varn = depvars[1].split('=')[1]
1934        qtot = ncobj.variables[depvars[2]]
1935        qtotv = ncobj.variables[depvars[2]]
1936        Nspecies = len(depvars) - 2
1937        for iv in range(Nspecies):
1938            if ncobj.variables.has_key(depvars[iv+2]):
1939                var1 = ncobj.variables[depvars[iv+2]]
1940                qtotv = qtotv + var1[:]
1941
1942        dnamesvar = list(qtot.dimensions)
1943        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1944
1945        diagout, diagoutd, diagoutvd = diag.compute_clwvl(var0, qtotv, dnamesvar,dvnamesvar)
1946
1947        # Removing the nonChecking variable-dimensions from the initial list
1948        varsadd = []
1949        diagoutvd = list(dvnames)
1950        for nonvd in NONchkvardims:
1951            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1952            varsadd.append(nonvd)
1953        ncvar.insert_variable(ncobj, varn, diagout, diagoutd, diagoutvd, newnc)
1954
1955# WRFgeop geopotential from WRF as PH + PHB
1956    elif diagn == 'WRFgeop':
1957        var0 = ncobj.variables[depvars[0]][:]
1958        var1 = ncobj.variables[depvars[1]][:]
1959
1960        # de-staggering geopotential
1961        diagout0 = var0 + var1
1962        dt = diagout0.shape[0]
1963        dz = diagout0.shape[1]
1964        dy = diagout0.shape[2]
1965        dx = diagout0.shape[3]
1966
1967        diagout = np.zeros((dt,dz-1,dy,dx), dtype=np.float)
1968        diagout = 0.5*(diagout0[:,1:dz,:,:]+diagout0[:,0:dz-1,:,:])
1969
1970        # Removing the nonChecking variable-dimensions from the initial list
1971        varsadd = []
1972        diagoutvd = list(dvnames)
1973        for nonvd in NONchkvardims:
1974            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1975            varsadd.append(nonvd)
1976
1977        ncvar.insert_variable(ncobj, 'zg', diagout, dnames, diagoutvd, newnc)
1978
1979# WRFpotevap_orPM potential evapotranspiration following Penman-Monteith formulation
1980#   implemented in ORCHIDEE (in src_sechiba/enerbil.f90) as: WRFdens, UST, U10, V10, T2, PSFC, QVAPOR
1981    elif diagn == 'WRFpotevap_orPM':
1982        var0 = WRFdens[:,0,:,:]
1983        var1 = ncobj.variables[depvars[1]][:]
1984        var2 = ncobj.variables[depvars[2]][:]
1985        var3 = ncobj.variables[depvars[3]][:]
1986        var4 = ncobj.variables[depvars[4]][:]
1987        var5 = ncobj.variables[depvars[5]][:]
1988        var6 = ncobj.variables[depvars[6]][:,0,:,:]
1989
1990        dnamesvar = list(ncobj.variables[depvars[1]].dimensions)
1991        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1992
1993        diagout = np.zeros(var1.shape, dtype=np.float)
1994        diagout, diagoutd, diagoutvd = diag.Forcompute_potevap_orPM(var0, var1, var2,\
1995          var3, var4, var5, var6, dnamesvar, dvnamesvar)
1996
1997        # Removing the nonChecking variable-dimensions from the initial list
1998        varsadd = []
1999        for nonvd in NONchkvardims:
2000            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2001            varsadd.append(nonvd)
2002
2003        ncvar.insert_variable(ncobj, 'evspsblpot', diagout, diagoutd, diagoutvd, newnc)
2004
2005# WRFmslp_ptarget sea-level pressure following ECMWF method as PSFC, HGT, WRFt, WRFp, ZNU, ZNW
2006    elif diagn == 'WRFpsl_ecmwf':
2007        var0 = ncobj.variables[depvars[0]][:]
2008        var1 = ncobj.variables[depvars[1]][0,:,:]
2009        var2 = WRFt[:,0,:,:]
2010        var4 = WRFp[:,0,:,:]
2011        var5 = ncobj.variables[depvars[4]][0,:]
2012        var6 = ncobj.variables[depvars[5]][0,:]
2013
2014        # This is quite too appriximate!! passing pressure at half-levels to 2nd full
2015        #   level, using eta values at full (ZNW) and half (ZNU) mass levels
2016        var3 = WRFp[:,0,:,:] + (var6[1] - var5[0])*(WRFp[:,1,:,:] - WRFp[:,0,:,:])/  \
2017          (var5[1]-var5[0])
2018
2019        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
2020        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2021
2022        diagout = np.zeros(var0.shape, dtype=np.float)
2023        diagout, diagoutd, diagoutvd = diag.Forcompute_psl_ecmwf(var0, var1, var2,   \
2024          var3, var4, dnamesvar, dvnamesvar)
2025
2026        # Removing the nonChecking variable-dimensions from the initial list
2027        varsadd = []
2028        for nonvd in NONchkvardims:
2029            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2030            varsadd.append(nonvd)
2031
2032        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
2033
2034# WRFmslp_ptarget sea-level pressure following ptarget method as WRFp, PSFC, WRFt, HGT, QVAPOR
2035    elif diagn == 'WRFpsl_ptarget':
2036        var0 = WRFp
2037        var1 = ncobj.variables[depvars[1]][:]
2038        var2 = WRFt
2039        var3 = ncobj.variables[depvars[3]][0,:,:]
2040        var4 = ncobj.variables[depvars[4]][:]
2041
2042        dnamesvar = list(ncobj.variables[depvars[4]].dimensions)
2043        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2044
2045        diagout = np.zeros(var0.shape, dtype=np.float)
2046        diagout, diagoutd, diagoutvd = diag.Forcompute_psl_ptarget(var0, var1, var2, \
2047          var3, var4, 700000., dnamesvar, dvnamesvar)
2048
2049        # Removing the nonChecking variable-dimensions from the initial list
2050        varsadd = []
2051        for nonvd in NONchkvardims:
2052            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2053            varsadd.append(nonvd)
2054
2055        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
2056
2057# WRFp pressure from WRF as P + PB
2058    elif diagn == 'WRFp':
2059        var0 = ncobj.variables[depvars[0]][:]
2060        var1 = ncobj.variables[depvars[1]][:]
2061           
2062        diagout = var0 + var1
2063        diagoutd = list(ncobj.variables[depvars[0]].dimensions)
2064        diagoutvd = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2065
2066        # Removing the nonChecking variable-dimensions from the initial list
2067        varsadd = []
2068        for nonvd in NONchkvardims:
2069            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2070            varsadd.append(nonvd)
2071
2072        ncvar.insert_variable(ncobj, 'pres', diagout, diagoutd, diagoutvd, newnc)
2073
2074# WRFpos
2075    elif diagn == 'WRFpos':
2076           
2077        dnamesvar = ncobj.variables['MAPFAC_M'].dimensions
2078        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2079
2080        ncvar.insert_variable(ncobj, 'WRFpos', WRFpos, dnamesvar, dvnamesvar, newnc)
2081
2082# WRFprw WRF water vapour path WRFdens, QVAPOR
2083    elif diagn == 'WRFprw':
2084           
2085        var0 = WRFdens
2086        var1 = ncobj.variables[depvars[1]]
2087
2088        dnamesvar = list(var1.dimensions)
2089        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2090
2091        diagout, diagoutd, diagoutvd = diag.compute_prw(var0, var1, dnamesvar,dvnamesvar)
2092
2093        # Removing the nonChecking variable-dimensions from the initial list
2094        varsadd = []
2095        diagoutvd = list(dvnames)
2096        for nonvd in NONchkvardims:
2097            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2098            varsadd.append(nonvd)
2099        ncvar.insert_variable(ncobj, 'prw', diagout, diagoutd, diagoutvd, newnc)
2100
2101# WRFrh (P, T, QVAPOR)
2102    elif diagn == 'WRFrh':
2103           
2104        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
2105        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2106
2107        ncvar.insert_variable(ncobj, 'hur', WRFrh, dnames, dvnames, newnc)
2108
2109# WRFrhs (PSFC, T2, Q2)
2110    elif diagn == 'WRFrhs':
2111           
2112        var0 = ncobj.variables[depvars[0]][:]
2113        var1 = ncobj.variables[depvars[1]][:]
2114        var2 = ncobj.variables[depvars[2]][:]
2115
2116        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
2117        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2118
2119        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
2120        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
2121
2122# rvors (u10, v10, WRFpos)
2123    elif diagn == 'WRFrvors':
2124           
2125        var0 = ncobj.variables[depvars[0]]
2126        var1 = ncobj.variables[depvars[1]]
2127
2128        diagout = rotational_z(var0, var1, distx)
2129
2130        dnamesvar = ncobj.variables[depvars[0]].dimensions
2131        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2132
2133        ncvar.insert_variable(ncobj, 'rvors', diagout, dnamesvar, dvnamesvar, newnc)
2134
2135# WRFt (T, P, PB)
2136    elif diagn == 'WRFt':
2137        var0 = ncobj.variables[depvars[0]][:]
2138        var1 = ncobj.variables[depvars[1]][:]
2139        var2 = ncobj.variables[depvars[2]][:]
2140
2141        p0=100000.
2142        p=var1 + var2
2143
2144        WRFt = (var0 + 300.)*(p/p0)**(2./7.)
2145
2146        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
2147        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2148
2149        # Removing the nonChecking variable-dimensions from the initial list
2150        varsadd = []
2151        diagoutvd = list(dvnames)
2152        for nonvd in NONchkvardims:
2153            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2154            varsadd.append(nonvd)
2155
2156        ncvar.insert_variable(ncobj, 'ta', WRFt, dnames, diagoutvd, newnc)
2157
2158# WRFtda (WRFrh, WRFt)
2159    elif diagn == 'WRFtda':
2160        ARM2 = fdef.module_definitions.arm2
2161        ARM3 = fdef.module_definitions.arm3
2162
2163        gammatarh = np.log(WRFrh) + ARM2*(WRFt-273.15)/((WRFt-273.15)+ARM3)
2164        td = ARM3*gammatarh/(ARM2-gammatarh)
2165
2166        dnamesvar = list(ncobj.variables['T'].dimensions)
2167        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2168
2169        # Removing the nonChecking variable-dimensions from the initial list
2170        varsadd = []
2171        diagoutvd = list(dvnames)
2172        for nonvd in NONchkvardims:
2173            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2174            varsadd.append(nonvd)
2175
2176        ncvar.insert_variable(ncobj, 'tda', td, dnames, diagoutvd, newnc)
2177
2178# WRFtdas (PSFC, T2, Q2)
2179    elif diagn == 'WRFtdas':
2180        ARM2 = fdef.module_definitions.arm2
2181        ARM3 = fdef.module_definitions.arm3
2182
2183        var0 = ncobj.variables[depvars[0]][:]
2184        var1 = ncobj.variables[depvars[1]][:]
2185        var2 = ncobj.variables[depvars[2]][:]
2186
2187        dnamesvar = list(ncobj.variables[depvars[1]].dimensions)
2188        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2189
2190        rhs, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
2191
2192        gammatarhs = np.log(rhs) + ARM2*(var1-273.15)/((var1-273.15)+ARM3)
2193        tdas = ARM3*gammatarhs/(ARM2-gammatarhs) + 273.15
2194
2195        # Removing the nonChecking variable-dimensions from the initial list
2196        varsadd = []
2197        diagoutvd = list(dvnames)
2198        for nonvd in NONchkvardims:
2199            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2200            varsadd.append(nonvd)
2201
2202        ncvar.insert_variable(ncobj, 'tdas', tdas, dnames, diagoutvd, newnc)
2203
2204# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
2205    elif diagn == 'WRFua':
2206        var0 = ncobj.variables[depvars[0]][:]
2207        var1 = ncobj.variables[depvars[1]][:]
2208        var2 = ncobj.variables[depvars[2]][:]
2209        var3 = ncobj.variables[depvars[3]][:]
2210
2211        # un-staggering variables
2212        if len(var0.shape) == 4:
2213            unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
2214        elif len(var0.shape) == 3:
2215            # Asuming sunding point (dimt, dimz, dimstgx)
2216            unstgdims = [var0.shape[0], var0.shape[1]]
2217
2218        ua = np.zeros(tuple(unstgdims), dtype=np.float)
2219        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
2220        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
2221
2222        if len(var0.shape) == 4:
2223            unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
2224            unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
2225
2226            for iz in range(var0.shape[1]):
2227                ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
2228
2229            dnamesvar = ['Time','bottom_top','south_north','west_east']
2230
2231        elif len(var0.shape) == 3:
2232            unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1])
2233            unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1])
2234            for iz in range(var0.shape[1]):
2235                ua[:,iz] = unstgvar0[:,iz]*var3 - unstgvar1[:,iz]*var2
2236
2237            dnamesvar = ['Time','bottom_top']
2238
2239        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2240
2241        # Removing the nonChecking variable-dimensions from the initial list
2242        varsadd = []
2243        diagoutvd = list(dvnames)
2244        for nonvd in NONchkvardims:
2245            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2246            varsadd.append(nonvd)
2247
2248        ncvar.insert_variable(ncobj, 'ua', ua, dnames, diagoutvd, newnc)
2249
2250# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
2251    elif diagn == 'WRFva':
2252        var0 = ncobj.variables[depvars[0]][:]
2253        var1 = ncobj.variables[depvars[1]][:]
2254        var2 = ncobj.variables[depvars[2]][:]
2255        var3 = ncobj.variables[depvars[3]][:]
2256
2257        # un-staggering variables
2258        if len(var0.shape) == 4:
2259            unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
2260        elif len(var0.shape) == 3:
2261            # Asuming sunding point (dimt, dimz, dimstgx)
2262            unstgdims = [var0.shape[0], var0.shape[1]]
2263
2264        va = np.zeros(tuple(unstgdims), dtype=np.float)
2265        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
2266        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
2267
2268        if len(var0.shape) == 4:
2269            unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
2270            unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
2271
2272            for iz in range(var0.shape[1]):
2273                va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
2274
2275            dnamesvar = ['Time','bottom_top','south_north','west_east']
2276
2277        elif len(var0.shape) == 3:
2278            unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1])
2279            unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1])
2280            for iz in range(var0.shape[1]):
2281                va[:,iz] = unstgvar0[:,iz]*var2 + unstgvar1[:,iz]*var3
2282
2283            dnamesvar = ['Time','bottom_top']
2284
2285        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2286
2287        # Removing the nonChecking variable-dimensions from the initial list
2288        varsadd = []
2289        diagoutvd = list(dvnames)
2290        for nonvd in NONchkvardims:
2291            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2292            varsadd.append(nonvd)
2293        ncvar.insert_variable(ncobj, 'va', va, dnames, diagoutvd, newnc)
2294
2295
2296# WRFwd (U, V, SINALPHA, COSALPHA) to be rotated !!
2297    elif diagn == 'WRFwd':
2298        var0 = ncobj.variables[depvars[0]][:]
2299        var1 = ncobj.variables[depvars[1]][:]
2300        var2 = ncobj.variables[depvars[2]][:]
2301        var3 = ncobj.variables[depvars[3]][:]
2302
2303        # un-staggering variables
2304        if len(var0.shape) == 4:
2305            unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
2306        elif len(var0.shape) == 3:
2307            # Asuming sunding point (dimt, dimz, dimstgx)
2308            unstgdims = [var0.shape[0], var0.shape[1]]
2309
2310        ua = np.zeros(tuple(unstgdims), dtype=np.float)
2311        va = np.zeros(tuple(unstgdims), dtype=np.float)
2312        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
2313        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
2314
2315        if len(var0.shape) == 4:
2316            unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
2317            unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
2318
2319            for iz in range(var0.shape[1]):
2320                ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
2321                va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
2322
2323            dnamesvar = ['Time','bottom_top','south_north','west_east']
2324
2325        elif len(var0.shape) == 3:
2326            unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1])
2327            unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1])
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']
2333
2334        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2335
2336        wd, dnames, dvnames = diag.compute_wd(ua, va, dnamesvar, dvnamesvar)
2337
2338        # Removing the nonChecking variable-dimensions from the initial list
2339        varsadd = []
2340        diagoutvd = list(dvnames)
2341        for nonvd in NONchkvardims:
2342            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2343            varsadd.append(nonvd)
2344
2345        ncvar.insert_variable(ncobj, 'wd', wd, dnames, diagoutvd, newnc)
2346
2347# WRFtime
2348    elif diagn == 'WRFtime':
2349           
2350        diagout = WRFtime
2351
2352        dnamesvar = ['Time']
2353        dvnamesvar = ['Times']
2354
2355        ncvar.insert_variable(ncobj, 'time', diagout, dnamesvar, dvnamesvar, newnc)
2356
2357# ws (U, V)
2358    elif diagn == 'ws':
2359           
2360        # un-staggering variables
2361        if len(var0.shape) == 4:
2362            unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
2363        elif len(var0.shape) == 3:
2364            # Asuming sunding point (dimt, dimz, dimstgx)
2365            unstgdims = [var0.shape[0], var0.shape[1]]
2366
2367        ua = np.zeros(tuple(unstgdims), dtype=np.float)
2368        va = np.zeros(tuple(unstgdims), dtype=np.float)
2369        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
2370        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
2371
2372        if len(var0.shape) == 4:
2373            unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
2374            unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
2375
2376            for iz in range(var0.shape[1]):
2377                ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
2378                va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
2379
2380            dnamesvar = ['Time','bottom_top','south_north','west_east']
2381
2382        elif len(var0.shape) == 3:
2383            unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1])
2384            unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1])
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']
2390
2391        diagout = np.sqrt(unstgvar0*unstgvar0 + unstgvar1*unstgvar1)
2392
2393#        dnamesvar = ncobj.variables[depvars[0]].dimensions
2394        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2395
2396        # Removing the nonChecking variable-dimensions from the initial list
2397        varsadd = []
2398        diagoutvd = list(dvnamesvar)
2399        for nonvd in NONchkvardims:
2400            if gen.searchInlist(dvnamesvar,nonvd): diagoutvd.remove(nonvd)
2401            varsadd.append(nonvd)
2402        ncvar.insert_variable(ncobj, 'ws', diagout, dnamesvar, diagoutvd, newnc)
2403
2404# wss (u10, v10)
2405    elif diagn == 'wss':
2406           
2407        var0 = ncobj.variables[depvars[0]][:]
2408        var1 = ncobj.variables[depvars[1]][:]
2409
2410        diagout = np.sqrt(var0*var0 + var1*var1)
2411
2412        dnamesvar = ncobj.variables[depvars[0]].dimensions
2413        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2414
2415        ncvar.insert_variable(ncobj, 'wss', diagout, dnamesvar, dvnamesvar, newnc)
2416
2417# WRFheight height from WRF geopotential as WRFGeop/g
2418    elif diagn == 'WRFheight':
2419           
2420        diagout = WRFgeop/grav
2421
2422        # Removing the nonChecking variable-dimensions from the initial list
2423        varsadd = []
2424        diagoutvd = list(dvnames)
2425        for nonvd in NONchkvardims:
2426            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2427            varsadd.append(nonvd)
2428
2429        ncvar.insert_variable(ncobj, 'zhgt', diagout, dnames, diagoutvd, newnc)
2430
2431# WRFheightrel relative-height from WRF geopotential as WRFgeop(PH + PHB)/g-HGT 'WRFheightrel|PH@PHB@HGT
2432    elif diagn == 'WRFheightrel':
2433        var0 = ncobj.variables[depvars[0]][:]
2434        var1 = ncobj.variables[depvars[1]][:]
2435        var2 = ncobj.variables[depvars[2]][:]
2436
2437        dimz = var0.shape[1]
2438        diagout = np.zeros(tuple(var0.shape), dtype=np.float)
2439        for iz in range(dimz):
2440            diagout[:,iz,:,:] = (var0[:,iz,:,:]+ var1[:,iz,:,:])/grav - var2
2441
2442        # Removing the nonChecking variable-dimensions from the initial list
2443        varsadd = []
2444        diagoutvd = list(dvnames)
2445        for nonvd in NONchkvardims:
2446            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2447            varsadd.append(nonvd)
2448
2449        ncvar.insert_variable(ncobj, 'zhgtrel', diagout, dnames, diagoutvd, newnc)
2450
2451# WRFzmla_gen generic boundary layer hieght computation from WRF theta, QVAPOR, WRFgeop, HGT,
2452    elif diagn == 'WRFzmlagen':
2453        var0 = ncobj.variables[depvars[0]][:]+300.
2454        var1 = ncobj.variables[depvars[1]][:]
2455        dimz = var0.shape[1]
2456        var2 = WRFgeop[:,1:dimz+1,:,:]/9.8
2457        var3 = ncobj.variables[depvars[3]][0,:,:]
2458
2459        diagout, diagoutd, diagoutvd = diag.Forcompute_zmla_gen(var0,var1,var2,var3, \
2460          dnames,dvnames)
2461
2462        # Removing the nonChecking variable-dimensions from the initial list
2463        varsadd = []
2464        for nonvd in NONchkvardims:
2465            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2466            varsadd.append(nonvd)
2467
2468        ncvar.insert_variable(ncobj, 'zmla', diagout, diagoutd, diagoutvd, newnc)
2469
2470# WRFzwind wind extrapolation at a given height using power law computation from WRF
2471#   U, V, WRFz, U10, V10, SINALPHA, COSALPHA, z=[zval]
2472    elif diagn == 'WRFzwind':
2473        var0 = ncobj.variables[depvars[0]][:]
2474        var1 = ncobj.variables[depvars[1]][:]
2475        var2 = WRFz
2476        var3 = ncobj.variables[depvars[3]][:]
2477        var4 = ncobj.variables[depvars[4]][:]
2478        var5 = ncobj.variables[depvars[5]][0,:,:]
2479        var6 = ncobj.variables[depvars[6]][0,:,:]
2480        var7 = np.float(depvars[7].split('=')[1])
2481
2482        # un-staggering 3D winds
2483        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
2484        va = np.zeros(tuple(unstgdims), dtype=np.float)
2485        unvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
2486        unvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
2487        unvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
2488        unvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
2489
2490        diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwind(unvar0,      \
2491          unvar1, var2, var3, var4, var5, var6, var7, dnames, dvnames)
2492
2493        # Removing the nonChecking variable-dimensions from the initial list
2494        varsadd = []
2495        for nonvd in NONchkvardims:
2496            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2497            varsadd.append(nonvd)
2498
2499        ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
2500        ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)
2501
2502# WRFzwind wind extrapolation at a given hieght using logarithmic law computation
2503#   from WRF U, V, WRFz, U10, V10, SINALPHA, COSALPHA, z=[zval]
2504    elif diagn == 'WRFzwind_log':
2505        var0 = ncobj.variables[depvars[0]][:]
2506        var1 = ncobj.variables[depvars[1]][:]
2507        var2 = WRFz
2508        var3 = ncobj.variables[depvars[3]][:]
2509        var4 = ncobj.variables[depvars[4]][:]
2510        var5 = ncobj.variables[depvars[5]][0,:,:]
2511        var6 = ncobj.variables[depvars[6]][0,:,:]
2512        var7 = np.float(depvars[7].split('=')[1])
2513
2514        # un-staggering 3D winds
2515        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
2516        va = np.zeros(tuple(unstgdims), dtype=np.float)
2517        unvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
2518        unvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
2519        unvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
2520        unvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
2521
2522        diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwind_log(unvar0,  \
2523          unvar1, var2, var3, var4, var5, var6, var7, dnames, dvnames)
2524
2525        # Removing the nonChecking variable-dimensions from the initial list
2526        varsadd = []
2527        for nonvd in NONchkvardims:
2528            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2529            varsadd.append(nonvd)
2530
2531        ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
2532        ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)
2533
2534# WRFzwindMO wind extrapolation at a given height computation using Monin-Obukhow
2535#   theory from WRF UST, ZNT, RMOL, U10, V10, SINALPHA, COSALPHA, z=[zval]
2536#   NOTE: only useful for [zval] < 80. m
2537    elif diagn == 'WRFzwindMO':
2538        var0 = ncobj.variables[depvars[0]][:]
2539        var1 = ncobj.variables[depvars[1]][:]
2540        var2 = ncobj.variables[depvars[2]][:]
2541        var3 = ncobj.variables[depvars[3]][:]
2542        var4 = ncobj.variables[depvars[4]][:]
2543        var5 = ncobj.variables[depvars[5]][0,:,:]
2544        var6 = ncobj.variables[depvars[6]][0,:,:]
2545        var7 = np.float(depvars[7].split('=')[1])
2546
2547        diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwindMO(var0, var1,\
2548          var2, var3, var4, var5, var6, var7, dnames, dvnames)
2549
2550        # Removing the nonChecking variable-dimensions from the initial list
2551        varsadd = []
2552        for nonvd in NONchkvardims:
2553            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2554            varsadd.append(nonvd)
2555
2556        ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
2557        ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)
2558
2559# zmla_gen generic boundary layer hieght computation from generic 2D-file theta, qv, zg, orog,
2560    elif diagn == 'zmlagen':
2561        var0 = ncobj.variables[depvars[0]][:]
2562        var1 = ncobj.variables[depvars[1]][:]
2563        dimz = var0.shape[1]
2564        var2 = ncobj.variables[depvars[2]][:]
2565        var3 = ncobj.variables[depvars[3]][:]
2566
2567        diagout, diagoutd, diagoutvd = diag.Forcompute_zmla_gen(var0,var1,var2,var3, \
2568          dnames,dvnames)
2569
2570        # Removing the nonChecking variable-dimensions from the initial list
2571        varsadd = []
2572        for nonvd in NONchkvardims:
2573            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2574            varsadd.append(nonvd)
2575
2576        ncvar.insert_variable(ncobj, 'zmla', diagout, diagoutd, diagoutvd, newnc)
2577
2578# zmla_genUWsnd generic boundary layer hieght computation from UWyoming sounding file theta, qv, zg
2579    elif diagn == 'zmlagenUWsnd':
2580        var0 = ncobj.variables[depvars[0]][:]
2581        var1 = ncobj.variables[depvars[1]][:]
2582        dimz = var0.shape[1]
2583        var2 = ncobj.variables[depvars[2]][:]/9.8
2584        var3 = ncobj.getncattr('Station_elevation')
2585
2586        diagout, diagoutd, diagoutvd = diag.Forcompute_zmla_gen(var0,var1,var2,var3, \
2587          dnames,dvnames)
2588
2589        # No need to remove topography height
2590        diagout = diagout + var3
2591
2592        # Removing the nonChecking variable-dimensions from the initial list
2593        varsadd = []
2594        for nonvd in NONchkvardims:
2595            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2596            varsadd.append(nonvd)
2597
2598        ncvar.insert_variable(ncobj, 'zmla', diagout, diagoutd, diagoutvd, newnc)
2599
2600    else:
2601        print errormsg
2602        print '  ' + main + ": diagnostic '" + diagn + "' not ready!!!"
2603        print '    available diagnostics: ', availdiags
2604        quit(-1)
2605
2606    newnc.sync()
2607    # Adding that additional variables required to compute some diagnostics which
2608    #   where not in the original file
2609    print '  adding additional variables...'
2610    for vadd in varsadd:
2611        if not gen.searchInlist(newnc.variables.keys(),vadd) and                     \
2612          dictcompvars.has_key(vadd):
2613            attrs = dictcompvars[vadd]
2614            vvn = attrs['name']
2615            if not gen.searchInlist(newnc.variables.keys(), vvn):
2616                iidvn = dvnames.index(vadd)
2617                dnn = dnames[iidvn]
2618                if vadd == 'WRFtime':
2619                    dvarvals = WRFtime[:]
2620                newvar = newnc.createVariable(vvn, 'f8', (dnn))
2621                newvar[:] = dvarvals
2622                for attn in attrs.keys():
2623                    if attn != 'name':
2624                        attv = attrs[attn]
2625                        ncvar.set_attribute(newvar, attn, attv)
2626
2627#   end of diagnostics
2628
2629# Global attributes
2630##
2631ncvar.add_global_PyNCplot(newnc, main, None, '2.0')
2632
2633gorigattrs = ncobj.ncattrs()
2634for attr in gorigattrs:
2635    attrv = ncobj.getncattr(attr)
2636    atvar = ncvar.set_attribute(newnc, attr, attrv)
2637
2638ncobj.close()
2639newnc.close()
2640
2641print '\n' + main + ': successfull writting of diagnostics file "' + ofile + '" !!!'
Note: See TracBrowser for help on using the repository browser.