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

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

Adding:

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