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

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

Adding:

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