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

Last change on this file since 1930 was 1927, checked in by lfita, 7 years ago

Adding

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