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

Last change on this file since 2049 was 2033, checked in by lfita, 7 years ago

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