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

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

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