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

Last change on this file since 2074 was 2072, checked in by lfita, 6 years ago

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