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

Last change on this file since 1802 was 1795, checked in by lfita, 7 years ago

Adding:

`psl_ecmwf': sea-level pressure computation following ECMWF

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