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

Last change on this file since 1697 was 1675, checked in by lfita, 8 years ago

Advancing diagnostics by creation of an independent module with all the variables called 'diag_tools.py'
File size: 39.8 KB
RevLine 
[1675]1# Python script to comput diagnostics
[365]2# L. Fita, LMD. CNR, UPMC-Jussieu, Paris, France
3# File diagnostics.inf provides the combination of variables to get the desired diagnostic
[772]4#   To be used with module_ForDiagnostics.F90, module_ForDiagnosticsVars.F90, module_generic.F90
[1150]5#      foudre: f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90 >& run_f2py.log
6#      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
[1149]7
[413]8## e.g. # diagnostics.py -d 'Time@time,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
9## e.g. # diagnostics.py -f /home/lluis/PY/diagnostics.inf -d variable_combo -v WRFprc
[365]10
[1675]11# Available general pupose diagnostics (model independent) providing (varv1, varv2, ..., dimns, dimvns)
12# compute_accum: Function to compute the accumulation of a variable
13# compute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following
14#   newmicro.F90 from LMDZ compute_clt(cldfra, pres, dimns, dimvns)
15# compute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ
16# compute_clivi: Function to compute cloud-ice water path (clivi)
17# compute_clwvl: Function to compute condensed water path (clwvl)
18# compute_deaccum: Function to compute the deaccumulation of a variable
19# compute_mslp: Function to compute mslp: mean sea level pressure following p_interp.F90 from WRF
20# compute_OMEGAw: Function to transform OMEGA [Pas-1] to velocities [ms-1]
21# compute_prw: Function to compute water vapour path (prw)
22# compute_rh: Function to compute relative humidity following 'Tetens' equation (T,P) ...'
23# compute_td: Function to compute the dew point temperature
24# compute_turbulence: Function to compute the rubulence term of the Taylor's decomposition ...'
25# compute_wds: Function to compute the wind direction
26# compute_wss: Function to compute the wind speed
27# compute_WRFuava: Function to compute geographical rotated WRF 3D winds
28# compute_WRFuasvas: Fucntion to compute geographical rotated WRF 2-meter winds
29# derivate_centered: Function to compute the centered derivate of a given field
30# def Forcompute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ via Fortran subroutine
31# Forcompute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ via a Fortran module
32
33# Others just providing variable values
34# var_cllmh: Fcuntion to compute cllmh on a 1D column
35# var_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ using 1D vertical column values
36# var_mslp: Fcuntion to compute mean sea-level pressure
37# var_virtualTemp: This function returns virtual temperature in K,
38# var_WRFtime: Function to copmute CFtimes from WRFtime variable
39# rotational_z: z-component of the rotatinoal of horizontal vectorial field
40# turbulence_var: Function to compute the Taylor's decomposition turbulence term from a a given variable
41
[365]42from optparse import OptionParser
43import numpy as np
44from netCDF4 import Dataset as NetCDFFile
45import os
46import re
47import nc_var_tools as ncvar
[756]48import generic_tools as gen
[654]49import datetime as dtime
[1163]50import module_ForDiag as fdin
[1675]51import diag_tools as diag
[365]52
53main = 'diagnostics.py'
54errormsg = 'ERROR -- error -- ERROR -- error'
55warnmsg = 'WARNING -- warning -- WARNING -- warning'
56
[654]57# Constants
58grav = 9.81
59
[365]60
61####### ###### ##### #### ### ## #
62comboinf="\nIF -d 'variable_combo', provides information of the combination to obtain -v [varn] with the ASCII file with the combinations as -f [combofile]"
63
64parser = OptionParser()
65parser.add_option("-f", "--netCDF_file", dest="ncfile", help="file to use", metavar="FILE")
66parser.add_option("-d", "--dimensions", dest="dimns", 
[1351]67  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)" + comboinf, 
[365]68  metavar="LABELS")
69parser.add_option("-v", "--variables", dest="varns", 
70  help=" [varn1]|[var11]@[...[varN1]],[...,[varnM]|[var1M]@[...[varLM]]] ',' list of variables to compute [varnK] and its necessary ones [var1K]...[varPK]", metavar="VALUES")
71
72(opts, args) = parser.parse_args()
73
74#######    #######
75## MAIN
76    #######
[649]77availdiags = ['ACRAINTOT', 'accum', 'clt', 'cllmh', 'deaccum', 'LMDZrh', 'mslp',     \
[1351]78  'OMEGAw', 'RAINTOT',                                                               \
[1581]79  'rvors', 'td', 'turbulence', 'WRFclivi', 'WRFclwvl', 'WRFgeop', 'WRFp',            \
80  'WRFrvors', 'ws', 'wds', 'wss', 'WRFheight', 'WRFheightrel', 'WRFua', 'WRFva']
[365]81
[649]82methods = ['accum', 'deaccum']
83
[365]84# Variables not to check
[612]85NONcheckingvars = ['cllmh', 'deaccum', 'TSrhs', 'TStd', 'TSwds', 'TSwss', 'WRFbils', \
[1581]86  'WRFclivi', 'WRFclwvl', 'WRFdens', 'WRFgeop',                                      \
[612]87  'WRFp', 'WRFtd',                                                                   \
[365]88  'WRFpos', 'WRFprc', 'WRFprls', 'WRFrh', 'LMDZrh', 'LMDZrhs', 'WRFrhs', 'WRFrvors', \
[914]89  'WRFt', 'WRFtime', 'WRFua', 'WRFva', 'WRFwds', 'WRFwss', 'WRFheight']
[365]90
[1351]91NONchkvardims = ['WRFtime']
92
[365]93ofile = 'diagnostics.nc'
94
95dimns = opts.dimns
96varns = opts.varns
97
98# Special method. knowing variable combination
99##
100if opts.dimns == 'variable_combo':
101    print warnmsg
102    print '  ' + main + ': knowing variable combination !!!'
103    combination = variable_combo(opts.varns,opts.ncfile)
104    print '     COMBO: ' + combination
105    quit(-1)
106
107if not os.path.isfile(opts.ncfile):
108    print errormsg
109    print '   ' + main + ": file '" + opts.ncfile + "' does not exist !!"
110    quit(-1)
111
112ncobj = NetCDFFile(opts.ncfile, 'r')
113
[1351]114# Looking for specific variables that might be use in more than one diagnostic
115WRFgeop_compute = False
116WRFp_compute = False
117WRFt_compute = False
118WRFrh_compute = False
119WRFght_compute = False
120WRFdens_compute = False
121WRFpos_compute = False
122WRFtime_compute = False
123
[365]124# File creation
125newnc = NetCDFFile(ofile,'w')
126
127# dimensions
128dimvalues = dimns.split(',')
129dnames = []
130dvnames = []
131
132for dimval in dimvalues:
[1351]133    dn = dimval.split('@')[0]
134    dnv = dimval.split('@')[1]
135    dnames.append(dn)
136    dvnames.append(dnv)
137    # Is there any dimension-variable which should be computed?
138    if dnv == 'WRFgeop':WRFgeop_compute = True
139    if dnv == 'WRFp': WRFp_compute = True
140    if dnv == 'WRFt': WRFt_compute = True
141    if dnv == 'WRFrh': WRFrh_compute = True
142    if dnv == 'WRFght': WRFght_compute = True
143    if dnv == 'WRFdens': WRFdens_compute = True
144    if dnv == 'WRFpos': WRFpos_compute = True
145    if dnv == 'WRFtime': WRFtime_compute = True
[365]146
147# diagnostics to compute
148diags = varns.split(',')
149Ndiags = len(diags)
150
151for idiag in range(Ndiags):
152    if diags[idiag].split('|')[1].find('@') == -1:
153        depvars = diags[idiag].split('|')[1]
[654]154        if depvars == 'WRFgeop':WRFgeop_compute = True
[365]155        if depvars == 'WRFp': WRFp_compute = True
156        if depvars == 'WRFt': WRFt_compute = True
157        if depvars == 'WRFrh': WRFrh_compute = True
158        if depvars == 'WRFght': WRFght_compute = True
159        if depvars == 'WRFdens': WRFdens_compute = True
160        if depvars == 'WRFpos': WRFpos_compute = True
[654]161        if depvars == 'WRFtime': WRFtime_compute = True
[365]162    else:
163        depvars = diags[idiag].split('|')[1].split('@')
[756]164        if gen.searchInlist(depvars, 'WRFgeop'): WRFgeop_compute = True
165        if gen.searchInlist(depvars, 'WRFp'): WRFp_compute = True
166        if gen.searchInlist(depvars, 'WRFt'): WRFt_compute = True
167        if gen.searchInlist(depvars, 'WRFrh'): WRFrh_compute = True
168        if gen.searchInlist(depvars, 'WRFght'): WRFght_compute = True
169        if gen.searchInlist(depvars, 'WRFdens'): WRFdens_compute = True
170        if gen.searchInlist(depvars, 'WRFpos'): WRFpos_compute = True
171        if gen.searchInlist(depvars, 'WRFtime'): WRFtime_compute = True
[365]172
[1351]173# Dictionary with the new computed variables to be able to add them
174dictcompvars = {}
[654]175if WRFgeop_compute:
176    print '  ' + main + ': Retrieving geopotential value from WRF as PH + PHB'
177    dimv = ncobj.variables['PH'].shape
178    WRFgeop = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
179
[1351]180    # Attributes of the variable
[1412]181    Vvals = gen.variables_values('WRFgeop')
[1351]182    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
183      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
184
[365]185if WRFp_compute:
186    print '  ' + main + ': Retrieving pressure value from WRF as P + PB'
187    dimv = ncobj.variables['P'].shape
188    WRFp = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
189
[1351]190    # Attributes of the variable
191    Vvals = gen.variables_values('WRFp')
192    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
193      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
194
[365]195if WRFght_compute:
196    print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...' 
197    WRFght = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
198
[1351]199    # Attributes of the variable
200    Vvals = gen.variables_values('WRFght')
201    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
202      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
203
[365]204if WRFrh_compute:
205    print '    ' + main + ": computing relative humidity from WRF as 'Tetens'" +     \
206      ' equation (T,P) ...'
207    p0=100000.
208    p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
209    tk = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
210    qv = ncobj.variables['QVAPOR'][:]
211
212    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
213    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
214
215    WRFrh = qv/data2
216
[1351]217    # Attributes of the variable
218    Vvals = gen.variables_values('WRFrh')
219    dictcompvars['WRFrh'] = {'name': Vvals[0], 'standard_name': Vvals[1],            \
220      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
221
[365]222if WRFt_compute:
223    print '    ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
224    p0=100000.
225    p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
226
227    WRFt = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
228
[1351]229    # Attributes of the variable
230    Vvals = gen.variables_values('WRFt')
231    dictcompvars['WRFt'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
232      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
233
[365]234if WRFdens_compute:
235    print '    ' + main + ': computing air density from WRF as ((MU + MUB) * ' +     \
236      'DNW)/g ...'
237
238# Just we need in in absolute values: Size of the central grid cell
239##    dxval = ncobj.getncattr('DX')
240##    dyval = ncobj.getncattr('DY')
241##    mapfac = ncobj.variables['MAPFAC_M'][:]
242##    area = dxval*dyval*mapfac
243
244    mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
245    dnw = ncobj.variables['DNW'][:]
246
247    WRFdens = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]),        \
248      dtype=np.float)
249    levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)
250
251    for it in range(mu.shape[0]):
252        for iz in range(dnw.shape[1]):
253            levval.fill(np.abs(dnw[it,iz]))
254            WRFdens[it,iz,:,:] = levval
255            WRFdens[it,iz,:,:] = mu[it,:,:]*WRFdens[it,iz,:,:]/grav
256
[1351]257    # Attributes of the variable
258    Vvals = gen.variables_values('WRFdens')
259    dictcompvars['WRFdens'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
260      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
261
[365]262if WRFpos_compute:
263# WRF positions from the lowest-leftest corner of the matrix
264    print '    ' + main + ': computing position from MAPFAC_M as sqrt(DY*j**2 + ' +  \
265      'DX*x**2)*MAPFAC_M ...'
266
267    mapfac = ncobj.variables['MAPFAC_M'][:]
268
269    distx = np.float(ncobj.getncattr('DX'))
270    disty = np.float(ncobj.getncattr('DY'))
271
272    print 'distx:',distx,'disty:',disty
273
274    dx = mapfac.shape[2]
275    dy = mapfac.shape[1]
276    dt = mapfac.shape[0]
277
278    WRFpos = np.zeros((dt, dy, dx), dtype=np.float)
279
280    for i in range(1,dx):
281        WRFpos[0,0,i] = distx*i/mapfac[0,0,i]
282    for j in range(1,dy):
283        i=0
284        WRFpos[0,j,i] = WRFpos[0,j-1,i] + disty/mapfac[0,j,i]
285        for i in range(1,dx):
286#            WRFpos[0,j,i] = np.sqrt((disty*j)**2. + (distx*i)**2.)/mapfac[0,j,i]
287#            WRFpos[0,j,i] = np.sqrt((disty*j)**2. + (distx*i)**2.)
288             WRFpos[0,j,i] = WRFpos[0,j,i-1] + distx/mapfac[0,j,i]
289
290    for it in range(1,dt):
291        WRFpos[it,:,:] = WRFpos[0,:,:]
292
[654]293if WRFtime_compute:
294    print '    ' + main + ': computing time from WRF as CFtime(Times) ...'
295
296    refdate='19491201000000'
297    tunitsval='minutes'
298
299    timeobj = ncobj.variables['Times']
300    timewrfv = timeobj[:]
301
302    yrref=refdate[0:4]
303    monref=refdate[4:6]
304    dayref=refdate[6:8]
305    horref=refdate[8:10]
306    minref=refdate[10:12]
307    secref=refdate[12:14]
308
309    refdateS = yrref + '-' + monref + '-' + dayref + ' ' + horref + ':' + minref +   \
310      ':' + secref
311
312    dt = timeobj.shape[0]
313    WRFtime = np.zeros((dt), dtype=np.float)
314
315    for it in range(dt):
[865]316        wrfdates = gen.datetimeStr_conversion(timewrfv[it,:],'WRFdatetime', 'matYmdHMS')
317        WRFtime[it] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
[654]318
319    tunits = tunitsval + ' since ' + refdateS
320
[1351]321    # Attributes of the variable
322    dictcompvars['WRFtime'] = {'name': 'time', 'standard_name': 'time',              \
323      'long_name': 'time', 'units': tunits, 'calendar': 'gregorian'}
324
[365]325### ## #
326# Going for the diagnostics
327### ## #
328print '  ' + main + ' ...'
[1404]329varsadd = []
[365]330
331for idiag in range(Ndiags):
332    print '    diagnostic:',diags[idiag]
333    diag = diags[idiag].split('|')[0]
334    depvars = diags[idiag].split('|')[1].split('@')
335    if diags[idiag].split('|')[1].find('@') != -1:
336        depvars = diags[idiag].split('|')[1].split('@')
337        if depvars[0] == 'deaccum': diag='deaccum'
[649]338        if depvars[0] == 'accum': diag='accum'
[365]339        for depv in depvars:
340            if not ncobj.variables.has_key(depv) and not                             \
[756]341              gen.searchInlist(NONcheckingvars, depv) and                            \
[865]342              not gen.searchInlist(methods, depv) and not depvars[0] == 'deaccum'    \
343              and not depvars[0] == 'accum':
[365]344                print errormsg
345                print '  ' + main + ": file '" + opts.ncfile +                       \
346                  "' does not have variable '" + depv + "' !!"
347                quit(-1)
348    else:
349        depvars = diags[idiag].split('|')[1]
350        if not ncobj.variables.has_key(depvars) and not                              \
[756]351          gen.searchInlist(NONcheckingvars, depvars) and                             \
352          not gen.searchInlist(methods, depvars):
353            print errormsg
[365]354            print '  ' + main + ": file '" + opts.ncfile +                           \
355              "' does not have variable '" + depvars + "' !!"
356            quit(-1)
357
358    print "\n    Computing '" + diag + "' from: ", depvars, '...'
359
360# acraintot: accumulated total precipitation from WRF RAINC, RAINNC
361    if diag == 'ACRAINTOT':
362           
363        var0 = ncobj.variables[depvars[0]]
364        var1 = ncobj.variables[depvars[1]]
365        diagout = var0[:] + var1[:]
366
367        dnamesvar = var0.dimensions
368        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
369
[1647]370        # Removing the nonChecking variable-dimensions from the initial list
371        varsadd = []
372        for nonvd in NONchkvardims:
373            if gen.searchInlist(dvnamesvar,nonvd): dvnamesvar.remove(nonvd)
374            varsadd.append(nonvd)
375
[649]376        ncvar.insert_variable(ncobj, 'pracc', diagout, dnamesvar, dvnamesvar, newnc)
[365]377
[649]378# accum: acumulation of any variable as (Variable, time [as [tunits]
379#   from/since ....], newvarname)
380    elif diag == 'accum':
381
382        var0 = ncobj.variables[depvars[0]]
383        var1 = ncobj.variables[depvars[1]]
384
385        dnamesvar = var0.dimensions
386        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
387
[1675]388        diagout, diagoutd, diagoutvd = diag.compute_accum(var0,dnamesvar,dvnamesvar)
[649]389
390        CFvarn = ncvar.variables_values(depvars[0])[0]
391
392# Removing the flux
393        if depvars[1] == 'XTIME':
394            dtimeunits = var1.getncattr('description')
395            tunits = dtimeunits.split(' ')[0]
396        else:
397            dtimeunits = var1.getncattr('units')
398            tunits = dtimeunits.split(' ')[0]
399
400        dtime = (var1[1] - var1[0])*timeunits_seconds(tunits)
401
402        ncvar.insert_variable(ncobj, CFvarn + 'acc', diagout*dtime, diagoutd, diagoutvd, newnc)
403
[365]404# cllmh with cldfra, pres
405    elif diag == 'cllmh':
406           
407        var0 = ncobj.variables[depvars[0]]
408        if depvars[1] == 'WRFp':
409            var1 = WRFp
410        else:
411            var01 = ncobj.variables[depvars[1]]
412            if len(size(var1.shape)) < len(size(var0.shape)):
413                var1 = np.brodcast_arrays(var01,var0)[0]
414            else:
415                var1 = var01
416
[1675]417        diagout, diagoutd, diagoutvd = diag.Forcompute_cllmh(var0,var1,dnames,dvnames)
[772]418
[1351]419        # Removing the nonChecking variable-dimensions from the initial list
420        varsadd = []
421        for nonvd in NONchkvardims:
422            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
423            varsadd.append(nonvd)
424
[365]425        ncvar.insert_variable(ncobj, 'cll', diagout[0,:], diagoutd, diagoutvd, newnc)
426        ncvar.insert_variable(ncobj, 'clm', diagout[1,:], diagoutd, diagoutvd, newnc)
427        ncvar.insert_variable(ncobj, 'clh', diagout[2,:], diagoutd, diagoutvd, newnc)
428
429# clt with cldfra
430    elif diag == 'clt':
431           
432        var0 = ncobj.variables[depvars]
[1675]433        diagout, diagoutd, diagoutvd = diag.Forcompute_clt(var0,dnames,dvnames)
[1351]434
435        # Removing the nonChecking variable-dimensions from the initial list
436        varsadd = []
437        for nonvd in NONchkvardims:
438            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
439            varsadd.append(nonvd)
440           
[365]441        ncvar.insert_variable(ncobj, 'clt', diagout, diagoutd, diagoutvd, newnc)
442
443# deaccum: deacumulation of any variable as (Variable, time [as [tunits]
444#   from/since ....], newvarname)
445    elif diag == 'deaccum':
446
447        var0 = ncobj.variables[depvars[1]]
448        var1 = ncobj.variables[depvars[2]]
449
450        dnamesvar = var0.dimensions
451        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
452
[1675]453        diagout, diagoutd, diagoutvd = diag.compute_deaccum(var0,dnamesvar,dvnamesvar)
[365]454
455# Transforming to a flux
456        if depvars[2] == 'XTIME':
457            dtimeunits = var1.getncattr('description')
458            tunits = dtimeunits.split(' ')[0]
459        else:
460            dtimeunits = var1.getncattr('units')
461            tunits = dtimeunits.split(' ')[0]
462
463        dtime = (var1[1] - var1[0])*timeunits_seconds(tunits)
464        ncvar.insert_variable(ncobj, depvars[3], diagout/dtime, diagoutd, diagoutvd, newnc)
465
466# LMDZrh (pres, t, r)
467    elif diag == 'LMDZrh':
468           
469        var0 = ncobj.variables[depvars[0]][:]
470        var1 = ncobj.variables[depvars[1]][:]
471        var2 = ncobj.variables[depvars[2]][:]
472
[1675]473        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnames,dvnames)
[1079]474        ncvar.insert_variable(ncobj, 'hur', diagout, diagoutd, diagoutvd, newnc)
[365]475
476# LMDZrhs (psol, t2m, q2m)
477    elif diag == 'LMDZrhs':
478           
479        var0 = ncobj.variables[depvars[0]][:]
480        var1 = ncobj.variables[depvars[1]][:]
481        var2 = ncobj.variables[depvars[2]][:]
482
483        dnamesvar = ncobj.variables[depvars[0]].dimensions
484        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
485
[1675]486        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
[365]487
[1079]488        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
[365]489
490# mslp: mean sea level pressure (pres, psfc, terrain, temp, qv)
491    elif diag == 'mslp' or diag == 'WRFmslp':
492           
493        var1 = ncobj.variables[depvars[1]][:]
494        var2 = ncobj.variables[depvars[2]][:]
495        var4 = ncobj.variables[depvars[4]][:]
496
497        if diag == 'WRFmslp':
498            var0 = WRFp
499            var3 = WRFt
500            dnamesvar = ncobj.variables['P'].dimensions
501        else:
502            var0 = ncobj.variables[depvars[0]][:]
503            var3 = ncobj.variables[depvars[3]][:]
504            dnamesvar = ncobj.variables[depvars[0]].dimensions
505
506        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
507
[1675]508        diagout, diagoutd, diagoutvd = diag.compute_mslp(var0, var1, var2, var3, var4,    \
[365]509          dnamesvar, dvnamesvar)
510
[1581]511        # Removing the nonChecking variable-dimensions from the initial list
512        varsadd = []
513        diagoutvd = list(dvnames)
514        for nonvd in NONchkvardims:
515            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
516            varsadd.append(nonvd)
[365]517        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
518
[642]519# OMEGAw (omega, p, t) from NCL formulation (https://www.ncl.ucar.edu/Document/Functions/Contributed/omega_to_w.shtml)
520    elif diag == 'OMEGAw':
521           
522        var0 = ncobj.variables[depvars[0]][:]
523        var1 = ncobj.variables[depvars[1]][:]
[643]524        var2 = ncobj.variables[depvars[2]][:]
[642]525
526        dnamesvar = ncobj.variables[depvars[0]].dimensions
527        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
528
[1675]529        diagout, diagoutd, diagoutvd = diag.compute_OMEGAw(var0,var1,var2,dnamesvar,dvnamesvar)
[642]530
531        ncvar.insert_variable(ncobj, 'wa', diagout, diagoutd, diagoutvd, newnc)
532
[365]533# raintot: instantaneous total precipitation from WRF as (RAINC + RAINC) / dTime
534    elif diag == 'RAINTOT':
535
536        var0 = ncobj.variables[depvars[0]]
537        var1 = ncobj.variables[depvars[1]]
[445]538        if depvars[2] != 'WRFtime':
[443]539            var2 = ncobj.variables[depvars[2]]
[654]540        else:
541            var2 = np.arange(var0.shape[0], dtype=int)
[365]542
543        var = var0[:] + var1[:]
544
545        dnamesvar = var0.dimensions
546        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
547
[1675]548        diagout, diagoutd, diagoutvd = diag.compute_deaccum(var,dnamesvar,dvnamesvar)
[365]549
550# Transforming to a flux
[654]551        if var2.shape[0] > 1:
[600]552            if depvars[2] != 'WRFtime':
553                dtimeunits = var2.getncattr('units')
554                tunits = dtimeunits.split(' ')[0]
555   
556                dtime = (var2[1] - var2[0])*timeunits_seconds(tunits)
557            else:
558                var2 = ncobj.variables['Times']
559                time1 = var2[0,:]
560                time2 = var2[1,:]
561                tmf1 = ''
562                tmf2 = ''
563                for ic in range(len(time1)):
564                    tmf1 = tmf1 + time1[ic]
565                    tmf2 = tmf2 + time2[ic]
[654]566                dtdate1 = dtime.datetime.strptime(tmf1,"%Y-%m-%d_%H:%M:%S")
567                dtdate2 = dtime.datetime.strptime(tmf2,"%Y-%m-%d_%H:%M:%S")
[600]568                diffdate12 = dtdate2 - dtdate1
569                dtime = diffdate12.total_seconds()
570                print 'dtime:',dtime
[442]571        else:
[600]572            print warnmsg
[1645]573            print '  ' + main + ": only 1 time-step for '" + diag + "' !!"
[600]574            print '    leaving a zero value!'
[1646]575            diagout = var0[:]*0.
[600]576            dtime=1.
[442]577
[1644]578        # Removing the nonChecking variable-dimensions from the initial list
579        varsadd = []
580        for nonvd in NONchkvardims:
581            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
582            varsadd.append(nonvd)
583           
[365]584        ncvar.insert_variable(ncobj, 'pr', diagout/dtime, diagoutd, diagoutvd, newnc)
585
[612]586# rhs (psfc, t, q) from TimeSeries files
587    elif diag == 'TSrhs':
588           
589        p0=100000.
590        var0 = ncobj.variables[depvars[0]][:]
591        var1 = (ncobj.variables[depvars[1]][:])*(var0/p0)**(2./7.)
592        var2 = ncobj.variables[depvars[2]][:]
593
594        dnamesvar = ncobj.variables[depvars[0]].dimensions
595        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
596
[1675]597        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
[612]598
[1079]599        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
[612]600
601# td (psfc, t, q) from TimeSeries files
[613]602    elif diag == 'TStd' or diag == 'td':
[612]603           
604        var0 = ncobj.variables[depvars[0]][:]
605        var1 = ncobj.variables[depvars[1]][:] - 273.15
606        var2 = ncobj.variables[depvars[2]][:]
607
608        dnamesvar = ncobj.variables[depvars[0]].dimensions
609        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
610
[1675]611        diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
[612]612
613        ncvar.insert_variable(ncobj, 'tds', diagout, diagoutd, diagoutvd, newnc)
614
615# td (psfc, t, q) from TimeSeries files
[616]616    elif diag == 'TStdC' or diag == 'tdC':
[612]617           
618        var0 = ncobj.variables[depvars[0]][:]
619# Temperature is already in degrees Celsius
620        var1 = ncobj.variables[depvars[1]][:]
621        var2 = ncobj.variables[depvars[2]][:]
622
623        dnamesvar = ncobj.variables[depvars[0]].dimensions
624        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
625
[1675]626        diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
[612]627
628        ncvar.insert_variable(ncobj, 'tds', diagout, diagoutd, diagoutvd, newnc)
629
630# wds (u, v)
631    elif diag == 'TSwds' or diag == 'wds' :
632 
633        var0 = ncobj.variables[depvars[0]][:]
634        var1 = ncobj.variables[depvars[1]][:]
635
636        dnamesvar = ncobj.variables[depvars[0]].dimensions
637        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
638
[1675]639        diagout, diagoutd, diagoutvd = diag.compute_wds(var0,var1,dnamesvar,dvnamesvar)
[612]640
641        ncvar.insert_variable(ncobj, 'wds', diagout, diagoutd, diagoutvd, newnc)
642
643# wss (u, v)
[613]644    elif diag == 'TSwss' or diag == 'wss':
[612]645           
646        var0 = ncobj.variables[depvars[0]][:]
647        var1 = ncobj.variables[depvars[1]][:]
648
649        dnamesvar = ncobj.variables[depvars[0]].dimensions
650        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
651
[1675]652        diagout, diagoutd, diagoutvd = diag.compute_wss(var0,var1,dnamesvar,dvnamesvar)
[612]653
654        ncvar.insert_variable(ncobj, 'wss', diagout, diagoutd, diagoutvd, newnc)
655
[365]656# turbulence (var)
657    elif diag == 'turbulence':
658
659        var0 = ncobj.variables[depvars][:]
660
661        dnamesvar = list(ncobj.variables[depvars].dimensions)
662        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
663
[1675]664        diagout, diagoutd, diagoutvd = diag.compute_turbulence(var0,dnamesvar,dvnamesvar)
[959]665        valsvar = gen.variables_values(depvars)
[365]666
[959]667        newvarn = depvars + 'turb'
668        print main + '; Lluis newvarn:', newvarn
669        ncvar.insert_variable(ncobj, newvarn, diagout, diagoutd, 
[365]670          diagoutvd, newnc)
[959]671        print main + '; Lluis variables:', newnc.variables.keys()
672        varobj = newnc.variables[newvarn]
[365]673        attrv = varobj.long_name
674        attr = varobj.delncattr('long_name')
675        newattr = ncvar.set_attribute(varobj, 'long_name', attrv +                   \
676          " Taylor decomposition turbulence term")
677
[390]678# WRFbils fom WRF as HFX + LH
679    elif diag == 'WRFbils':
680           
681        var0 = ncobj.variables[depvars[0]][:]
682        var1 = ncobj.variables[depvars[1]][:]
683
684        diagout = var0 + var1
[867]685        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
686        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
[390]687
[867]688        ncvar.insert_variable(ncobj, 'bils', diagout, dnamesvar, dvnamesvar, newnc)
[390]689
[1581]690# WRFclivi WRF water vapour path WRFdens, QICE, QGRAUPEL, QHAIL
691    elif diag == 'WRFclivi':
692           
693        var0 = WRFdens
694        qtot = ncobj.variables[depvars[1]]
695        qtotv = qtot[:]
696        Nspecies = len(depvars) - 2
697        for iv in range(Nspecies):
[1585]698            if ncobj.variables.has_key(depvars[iv+2]):
699                var1 = ncobj.variables[depvars[iv+2]][:]
700                qtotv = qtotv + var1
[1581]701
702        dnamesvar = list(qtot.dimensions)
703        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
704
[1675]705        diagout, diagoutd, diagoutvd = diag.compute_clivi(var0, qtotv, dnamesvar,dvnamesvar)
[1581]706
707        # Removing the nonChecking variable-dimensions from the initial list
708        varsadd = []
709        diagoutvd = list(dvnames)
710        for nonvd in NONchkvardims:
711            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
712            varsadd.append(nonvd)
713        ncvar.insert_variable(ncobj, 'clivi', diagout, diagoutd, diagoutvd, newnc)
714
715# WRFclwvl WRF water cloud-condensed path WRFdens, QCLOUD, QICE, QGRAUPEL, QHAIL
716    elif diag == 'WRFclwvl':
717           
718        var0 = WRFdens
719        qtot = ncobj.variables[depvars[1]]
720        qtotv = ncobj.variables[depvars[1]]
721        Nspecies = len(depvars) - 2
722        for iv in range(Nspecies):
[1585]723            if ncobj.variables.has_key(depvars[iv+2]):
724                var1 = ncobj.variables[depvars[iv+2]]
725                qtotv = qtotv + var1[:]
[1581]726
727        dnamesvar = list(qtot.dimensions)
728        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
729
[1675]730        diagout, diagoutd, diagoutvd = diag.compute_clwvl(var0, qtotv, dnamesvar,dvnamesvar)
[1581]731
732        # Removing the nonChecking variable-dimensions from the initial list
733        varsadd = []
734        diagoutvd = list(dvnames)
735        for nonvd in NONchkvardims:
736            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
737            varsadd.append(nonvd)
738        ncvar.insert_variable(ncobj, 'clwvl', diagout, diagoutd, diagoutvd, newnc)
739
[654]740# WRFgeop geopotential from WRF as PH + PHB
741    elif diag == 'WRFgeop':
[1382]742        var0 = ncobj.variables[depvars[0]][:]
743        var1 = ncobj.variables[depvars[1]][:]
[654]744
[1382]745        # de-staggering geopotential
746        diagout0 = var0 + var1
747        dt = diagout0.shape[0]
748        dz = diagout0.shape[1]
749        dy = diagout0.shape[2]
750        dx = diagout0.shape[3]
751
752        diagout = np.zeros((dt,dz-1,dy,dx), dtype=np.float)
753        diagout = 0.5*(diagout0[:,1:dz,:,:]+diagout0[:,0:dz-1,:,:])
754
755        # Removing the nonChecking variable-dimensions from the initial list
756        varsadd = []
[1389]757        diagoutvd = list(dvnames)
[1382]758        for nonvd in NONchkvardims:
[1389]759            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
[1382]760            varsadd.append(nonvd)
761
[1389]762        ncvar.insert_variable(ncobj, 'zg', diagout, dnames, diagoutvd, newnc)
[654]763
[390]764# WRFp pressure from WRF as P + PB
[365]765    elif diag == 'WRFp':
766           
767        diagout = WRFp
768
769        ncvar.insert_variable(ncobj, 'pres', diagout, dnames, dvnames, newnc)
770
771# WRFpos
772    elif diag == 'WRFpos':
773           
774        dnamesvar = ncobj.variables['MAPFAC_M'].dimensions
775        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
776
777        ncvar.insert_variable(ncobj, 'WRFpos', WRFpos, dnamesvar, dvnamesvar, newnc)
778
779# WRFprw WRF water vapour path WRFdens, QVAPOR
780    elif diag == 'WRFprw':
781           
782        var0 = WRFdens
783        var1 = ncobj.variables[depvars[1]]
784
785        dnamesvar = list(var1.dimensions)
786        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
787
[1675]788        diagout, diagoutd, diagoutvd = diag.compute_prw(var0, var1, dnamesvar,dvnamesvar)
[365]789
[1586]790        # Removing the nonChecking variable-dimensions from the initial list
791        varsadd = []
792        diagoutvd = list(dvnames)
793        for nonvd in NONchkvardims:
794            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
795            varsadd.append(nonvd)
[365]796        ncvar.insert_variable(ncobj, 'prw', diagout, diagoutd, diagoutvd, newnc)
797
798# WRFrh (P, T, QVAPOR)
799    elif diag == 'WRFrh':
800           
801        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
802        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
803
[878]804        ncvar.insert_variable(ncobj, 'hur', WRFrh, dnames, dvnames, newnc)
[365]805
806# WRFrhs (PSFC, T2, Q2)
807    elif diag == 'WRFrhs':
808           
809        var0 = ncobj.variables[depvars[0]][:]
810        var1 = ncobj.variables[depvars[1]][:]
811        var2 = ncobj.variables[depvars[2]][:]
812
813        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
814        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
815
[1675]816        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
[878]817        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
[365]818
819# rvors (u10, v10, WRFpos)
820    elif diag == 'WRFrvors':
821           
822        var0 = ncobj.variables[depvars[0]]
823        var1 = ncobj.variables[depvars[1]]
824
825        diagout = rotational_z(var0, var1, distx)
826
827        dnamesvar = ncobj.variables[depvars[0]].dimensions
828        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
829
830        ncvar.insert_variable(ncobj, 'rvors', diagout, dnamesvar, dvnamesvar, newnc)
831
[884]832# WRFt (T, P, PB)
833    elif diag == 'WRFt':
834        var0 = ncobj.variables[depvars[0]][:]
835        var1 = ncobj.variables[depvars[1]][:]
836        var2 = ncobj.variables[depvars[2]][:]
[654]837
[884]838        p0=100000.
839        p=var1 + var2
840
841        WRFt = (var0 + 300.)*(p/p0)**(2./7.)
842
843        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
844        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
845
[1382]846        # Removing the nonChecking variable-dimensions from the initial list
847        varsadd = []
[1389]848        diagoutvd = list(dvnames)
[1382]849        for nonvd in NONchkvardims:
[1389]850            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
[1382]851            varsadd.append(nonvd)
852
[1389]853        ncvar.insert_variable(ncobj, 'ta', WRFt, dnames, diagoutvd, newnc)
[884]854
[914]855# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
856    elif diag == 'WRFua':
857        var0 = ncobj.variables[depvars[0]][:]
858        var1 = ncobj.variables[depvars[1]][:]
859        var2 = ncobj.variables[depvars[2]][:]
860        var3 = ncobj.variables[depvars[3]][:]
861
862        # un-staggering variables
863        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
864        ua = np.zeros(tuple(unstgdims), dtype=np.float)
865        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
866        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
867        unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
868        unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
869
870        for iz in range(var0.shape[1]):
871            ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
872
[959]873        dnamesvar = ['Time','bottom_top','south_north','west_east']
[914]874        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
875
[1404]876        # Removing the nonChecking variable-dimensions from the initial list
877        varsadd = []
878        diagoutvd = list(dvnames)
879        for nonvd in NONchkvardims:
880            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
881            varsadd.append(nonvd)
[914]882
[1404]883        ncvar.insert_variable(ncobj, 'ua', ua, dnames, diagoutvd, newnc)
884
[914]885# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
886    elif diag == 'WRFva':
887        var0 = ncobj.variables[depvars[0]][:]
888        var1 = ncobj.variables[depvars[1]][:]
889        var2 = ncobj.variables[depvars[2]][:]
890        var3 = ncobj.variables[depvars[3]][:]
891
892        # un-staggering variables
893        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
894        va = np.zeros(tuple(unstgdims), dtype=np.float)
895        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
896        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
897        unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
898        unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
899        for iz in range(var0.shape[1]):
900            va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
901
[959]902        dnamesvar = ['Time','bottom_top','south_north','west_east']
[914]903        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
904
[1404]905        # Removing the nonChecking variable-dimensions from the initial list
906        varsadd = []
907        diagoutvd = list(dvnames)
908        for nonvd in NONchkvardims:
909            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
910            varsadd.append(nonvd)
911        ncvar.insert_variable(ncobj, 'va', va, dnames, diagoutvd, newnc)
[914]912
913# WRFtime
[654]914    elif diag == 'WRFtime':
915           
916        diagout = WRFtime
917
918        dnamesvar = ['Time']
919        dvnamesvar = ['Times']
920
921        ncvar.insert_variable(ncobj, 'time', diagout, dnamesvar, dvnamesvar, newnc)
922
[959]923# ws (U, V)
924    elif diag == 'ws':
925           
926        var0 = ncobj.variables[depvars[0]][:]
927        var1 = ncobj.variables[depvars[1]][:]
928        # un-staggering variables
929        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
930        va = np.zeros(tuple(unstgdims), dtype=np.float)
931        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
932        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
933        unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
934        unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
935
936        dnamesvar = ['Time','bottom_top','south_north','west_east']
937        diagout = np.sqrt(unstgvar0*unstgvar0 + unstgvar1*unstgvar1)
938
939#        dnamesvar = ncobj.variables[depvars[0]].dimensions
940        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
941
[1408]942        # Removing the nonChecking variable-dimensions from the initial list
943        varsadd = []
944        diagoutvd = list(dvnamesvar)
945        for nonvd in NONchkvardims:
946            if gen.searchInlist(dvnamesvar,nonvd): diagoutvd.remove(nonvd)
947            varsadd.append(nonvd)
948        ncvar.insert_variable(ncobj, 'ws', diagout, dnamesvar, diagoutvd, newnc)
[959]949
[365]950# wss (u10, v10)
951    elif diag == 'wss':
952           
953        var0 = ncobj.variables[depvars[0]][:]
954        var1 = ncobj.variables[depvars[1]][:]
955
956        diagout = np.sqrt(var0*var0 + var1*var1)
957
958        dnamesvar = ncobj.variables[depvars[0]].dimensions
959        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
960
961        ncvar.insert_variable(ncobj, 'wss', diagout, dnamesvar, dvnamesvar, newnc)
962
[654]963# WRFheight height from WRF geopotential as WRFGeop/g
964    elif diag == 'WRFheight':
965           
966        diagout = WRFgeop/grav
967
[1412]968        # Removing the nonChecking variable-dimensions from the initial list
969        varsadd = []
970        diagoutvd = list(dvnames)
971        for nonvd in NONchkvardims:
972            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
973            varsadd.append(nonvd)
[654]974
[1412]975        ncvar.insert_variable(ncobj, 'zhgt', diagout, dnames, diagoutvd, newnc)
976
[1413]977# WRFheightrel relative-height from WRF geopotential as WRFgeop(PH + PHB)/g-HGT 'WRFheightrel|PH@PHB@HGT
978    elif diag == 'WRFheightrel':
979        var0 = ncobj.variables[depvars[0]][:]
980        var1 = ncobj.variables[depvars[1]][:]
981        var2 = ncobj.variables[depvars[2]][:]
982
983        dimz = var0.shape[1]
984        diagout = np.zeros(tuple(var0.shape), dtype=np.float)
985        for iz in range(dimz):
[1419]986            diagout[:,iz,:,:] = (var0[:,iz,:,:]+ var1[:,iz,:,:])/grav - var2
[1413]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
995        ncvar.insert_variable(ncobj, 'zhgtrel', diagout, dnames, diagoutvd, newnc)
996
[365]997    else:
998        print errormsg
999        print '  ' + main + ": diagnostic '" + diag + "' not ready!!!"
1000        print '    available diagnostics: ', availdiags
1001        quit(-1)
1002
1003    newnc.sync()
[1351]1004    # Adding that additional variables required to compute some diagnostics which
1005    #   where not in the original file
1006    for vadd in varsadd:
1007        if not gen.searchInlist(newnc.variables.keys(),vadd):
1008            attrs = dictcompvars[vadd]
1009            vvn = attrs['name']
1010            if not gen.searchInlist(newnc.variables.keys(), vvn):
1011                iidvn = dvnames.index(vadd)
1012                dnn = dnames[iidvn]
1013                if vadd == 'WRFtime':
1014                    dvarvals = WRFtime[:]
1015                newvar = newnc.createVariable(vvn, 'f8', (dnn))
1016                newvar[:] = dvarvals
1017                for attn in attrs.keys():
1018                    if attn != 'name':
1019                        attv = attrs[attn]
1020                        ncvar.set_attribute(newvar, attn, attv)
[365]1021
1022#   end of diagnostics
1023
1024# Global attributes
1025##
[1675]1026add_global_PyNCplot(newnc, main, None, '2.0')
[365]1027
1028gorigattrs = ncobj.ncattrs()
1029for attr in gorigattrs:
1030    attrv = ncobj.getncattr(attr)
1031    atvar = ncvar.set_attribute(newnc, attr, attrv)
1032
1033ncobj.close()
1034newnc.close()
1035
1036print '\n' + main + ': successfull writting of diagnostics file "' + ofile + '" !!!'
Note: See TracBrowser for help on using the repository browser.