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

Last change on this file since 2235 was 2223, checked in by lfita, 7 years ago

Adding to 'range_faces':

  • -2:|downhill;-1:|downhill-valley;|0:|peak|;1:|uphill-valley;2:|uphill
  • filtering outside any range
File size: 79.6 KB
Line 
1# Python script to comput diagnostics
2# From L. Fita work in different places: CCRC (Australia), LMD (France)
3# More information at: http://www.xn--llusfb-5va.cat/python/PyNCplot
4#
5# pyNCplot and its component nc_var.py comes with ABSOLUTELY NO WARRANTY.
6# This work is licendes under a Creative Commons
7#   Attribution-ShareAlike 4.0 International License (http://creativecommons.org/licenses/by-sa/4.0)
8#
9# L. Fita, CIMA. CONICET-UBA, CNRS UMI-IFAECI, C.A. Buenos Aires, Argentina
10# File diagnostics.inf provides the combination of variables to get the desired diagnostic
11#   To be used with module_ForDiagnostics.F90, module_ForDiagnosticsVars.F90, module_generic.F90
12#      foudre: f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90 >& run_f2py.log
13#      ciclad: f2py --f90flags="-fPIC" --f90exec=/usr/bin/gfortran -L/opt/canopy-1.3.0/Canopy_64bit/System/lib/ -L/usr/lib64/ -L/opt/canopy-1.3.0/Canopy_64bit/System/lib/ -m module_ForDiagnostics -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90 >& run_f2py.log
14
15## e.g. # diagnostics.py -d 'Time@WRFtime,bottom_top@ZNU,south_north@XLAT,west_east@XLONG' -v 'clt|CLDFRA,cllmh|CLDFRA@WRFp,RAINTOT|RAINC@RAINNC@RAINSH@XTIME' -f WRF_LMDZ/NPv31/wrfout_d01_1980-03-01_00:00:00
16## e.g. # diagnostics.py -f /home/lluis/PY/diagnostics.inf -d variable_combo -v WRFprc
17
18# Available general pupose diagnostics (model independent) providing (varv1, varv2, ..., dimns, dimvns)
19# compute_accum: Function to compute the accumulation of a variable
20# compute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following
21#   newmicro.F90 from LMDZ compute_clt(cldfra, pres, dimns, dimvns)
22# compute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ
23# compute_clivi: Function to compute cloud-ice water path (clivi)
24# compute_clwvl: Function to compute condensed water path (clwvl)
25# compute_deaccum: Function to compute the deaccumulation of a variable
26# compute_mslp: Function to compute mslp: mean sea level pressure following p_interp.F90 from WRF
27# compute_OMEGAw: Function to transform OMEGA [Pas-1] to velocities [ms-1]
28# compute_prw: Function to compute water vapour path (prw)
29# compute_range_faces: Function to compute faces [uphill, valley, downhill] of sections of a mountain
30#   range, along a given face
31# compute_rh: Function to compute relative humidity following 'Tetens' equation (T,P) ...'
32# compute_td: Function to compute the dew point temperature
33# compute_turbulence: Function to compute the rubulence term of the Taylor's decomposition ...'
34# compute_wds: Function to compute the wind direction
35# compute_wss: Function to compute the wind speed
36# compute_WRFuava: Function to compute geographical rotated WRF 3D winds
37# compute_WRFuasvas: Fucntion to compute geographical rotated WRF 2-meter winds
38# derivate_centered: Function to compute the centered derivate of a given field
39# def Forcompute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ via Fortran subroutine
40# Forcompute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ via a Fortran module
41# Forcompute_psl_ptarget: Function to compute the sea-level pressure following target_pressure value found in `p_interp.F'
42
43# Others just providing variable values
44# var_cllmh: Fcuntion to compute cllmh on a 1D column
45# var_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ using 1D vertical column values
46# var_mslp: Fcuntion to compute mean sea-level pressure
47# var_virtualTemp: This function returns virtual temperature in K,
48# var_WRFtime: Function to copmute CFtimes from WRFtime variable
49# rotational_z: z-component of the rotatinoal of horizontal vectorial field
50# turbulence_var: Function to compute the Taylor's decomposition turbulence term from a a given variable
51# timeoverthres: When a given variable [varname] overpass a given [value]. Being [CFvarn] the name of the diagnostics in
52#   variables_values.dat
53# timemax ([varname], time). When a given variable [varname] got its maximum
54
55from optparse import OptionParser
56import numpy as np
57from netCDF4 import Dataset as NetCDFFile
58import os
59import re
60import nc_var_tools as ncvar
61import generic_tools as gen
62import datetime as dtime
63import module_ForDiag as fdin
64import module_ForDef as fdef
65import diag_tools as diag
66
67main = 'diagnostics.py'
68errormsg = 'ERROR -- error -- ERROR -- error'
69warnmsg = 'WARNING -- warning -- WARNING -- warning'
70
71# Constants
72grav = 9.81
73
74
75####### ###### ##### #### ### ## #
76comboinf="\nIF -d 'variable_combo', provides information of the combination to obtain -v [varn] with the ASCII file with the combinations as -f [combofile]"
77
78parser = OptionParser()
79parser.add_option("-f", "--netCDF_file", dest="ncfile", help="file to use", metavar="FILE")
80parser.add_option("-d", "--dimensions", dest="dimns", 
81  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, 
82  metavar="LABELS")
83parser.add_option("-v", "--variables", dest="varns", 
84  help=" [varn1]|[var11]@[...[varN1]],[...,[varnM]|[var1M]@[...[varLM]]] ',' list of variables to compute [varnK] and its necessary ones [var1K]...[varPK]", metavar="VALUES")
85
86(opts, args) = parser.parse_args()
87
88#######    #######
89## MAIN
90    #######
91availdiags = ['ACRAINTOT', 'accum', 'clt', 'cllmh', 'convini', 'deaccum', 'fog_K84', \
92  'fog_RUC', 'LMDZrh', 'mslp', 'OMEGAw', 'RAINTOT', 'range_faces',                   \
93  'rvors', 'td', 'timemax', 'timeoverthres', 'turbulence', 'uavaFROMwswd',           \
94  'WRFcape_afwa', 'WRFclivi', 'WRFclwvi', 'WRF_denszint', 'WRFgeop',                 \
95  'WRFmrso', 'WRFmrsos', 'WRFpotevap_orPM', 'WRFp', 'WRFpsl_ecmwf',                  \
96  'WRFpsl_ptarget', 'WRFrvors', 'WRFslw', 'ws', 'wds', 'wss', 'WRFheight',           \
97  'WRFheightrel', 'WRFtda', 'WRFtdas', 'WRFua', 'WRFva', 'WRFzwind', 'WRFzwind_log', \
98  'WRFzwindMO']
99
100methods = ['accum', 'deaccum']
101
102# Variables not to check
103NONcheckingvars = ['accum', 'cllmh', 'deaccum', 'face', 'LONLATdxdy', 'TSrhs',       \
104  'TStd', 'TSwds', 'TSwss',                                                          \
105  'WRFbils',                                                                         \
106  'WRFclivi', 'WRFclwvi', 'WRFdens', 'WRFdx', 'WRFdxdy', 'WRFdy', 'WRFgeop',         \
107  'WRFp', 'WRFtd',                                                                   \
108  'WRFpos', 'WRFprc', 'WRFprls', 'WRFrh', 'LMDZrh', 'LMDZrhs',                       \
109  'WRFrhs', 'WRFrvors',                                                              \
110  'WRFt', 'WRFtime', 'WRFua', 'WRFva', 'WRFwds', 'WRFwss', 'WRFheight', 'WRFz',      \
111  'WRFzg']
112
113# diagnostics not to check their dependeny
114NONcheckdepvars = ['accum', 'deaccum', 'timeoverthres', 'WRF_denszint',              \
115  'WRFzwind_log', 'WRFzwind', 'WRFzwindMO']
116
117NONchkvardims = ['WRFtime']
118
119ofile = 'diagnostics.nc'
120
121dimns = opts.dimns
122varns = opts.varns
123
124# Special method. knowing variable combination
125##
126if opts.dimns == 'variable_combo':
127    print warnmsg
128    print '  ' + main + ': knowing variable combination !!!'
129    combination = variable_combo(opts.varns,opts.ncfile)
130    print '     COMBO: ' + combination
131    quit(-1)
132
133if opts.ncfile is None:
134    print errormsg
135    print '   ' + main + ": No file provided !!"
136    print '     is mandatory to provide a file -f [filename]'
137    quit(-1)
138
139if opts.dimns is None:
140    print errormsg
141    print '   ' + main + ": No description of dimensions are provided !!"
142    print '     is mandatory to provide description of dimensions as -d [dimn]@[vardimname],... '
143    quit(-1)
144
145if opts.varns is None:
146    print errormsg
147    print '   ' + main + ": No variable to diagnose is provided !!"
148    print '     is mandatory to provide a variable to diagnose as -v [diagn]|[varn1]@... '
149    quit(-1)
150
151if not os.path.isfile(opts.ncfile):
152    print errormsg
153    print '   ' + main + ": file '" + opts.ncfile + "' does not exist !!"
154    quit(-1)
155
156ncobj = NetCDFFile(opts.ncfile, 'r')
157
158# Looking for specific variables that might be use in more than one diagnostic
159WRFgeop_compute = False
160WRFp_compute = False
161WRFt_compute = False
162WRFrh_compute = False
163WRFght_compute = False
164WRFdens_compute = False
165WRFpos_compute = False
166WRFtime_compute = False
167WRFz_compute = False
168WRFdxdy_compute = False
169LONLATdxdy_compute = False
170
171# File creation
172newnc = NetCDFFile(ofile,'w')
173
174# dimensions
175dimvalues = dimns.split(',')
176dnames = []
177dvnames = []
178
179for dimval in dimvalues:
180    dn = dimval.split('@')[0]
181    dnv = dimval.split('@')[1]
182    dnames.append(dn)
183    dvnames.append(dnv)
184    # Is there any dimension-variable which should be computed?
185    if dnv == 'WRFgeop':WRFgeop_compute = True
186    if dnv == 'WRFp': WRFp_compute = True
187    if dnv == 'WRFt': WRFt_compute = True
188    if dnv == 'WRFrh': WRFrh_compute = True
189    if dnv == 'WRFght': WRFght_compute = True
190    if dnv == 'WRFdens': WRFdens_compute = True
191    if dnv == 'WRFpos': WRFpos_compute = True
192    if dnv == 'WRFtime': WRFtime_compute = True
193    if dnv == 'WRFz':WRFz_compute = True
194    if dnv == 'WRFdxdy':WRFdxdy_compute = True
195    if dnv == 'LONLATdxdy':LONLATdxdy_compute = True
196
197# diagnostics to compute
198diags = varns.split(',')
199Ndiags = len(diags)
200
201for idiag in range(Ndiags):
202    if diags[idiag].split('|')[1].find('@') == -1:
203        depvars = diags[idiag].split('|')[1]
204        if depvars == 'WRFgeop':WRFgeop_compute = True
205        if depvars == 'WRFp': WRFp_compute = True
206        if depvars == 'WRFt': WRFt_compute = True
207        if depvars == 'WRFrh': WRFrh_compute = True
208        if depvars == 'WRFght': WRFght_compute = True
209        if depvars == 'WRFdens': WRFdens_compute = True
210        if depvars == 'WRFpos': WRFpos_compute = True
211        if depvars == 'WRFtime': WRFtime_compute = True
212        if depvars == 'WRFz': WRFz_compute = True
213    else:
214        depvars = diags[idiag].split('|')[1].split('@')
215        if gen.searchInlist(depvars, 'WRFgeop'): WRFgeop_compute = True
216        if gen.searchInlist(depvars, 'WRFp'): WRFp_compute = True
217        if gen.searchInlist(depvars, 'WRFt'): WRFt_compute = True
218        if gen.searchInlist(depvars, 'WRFrh'): WRFrh_compute = True
219        if gen.searchInlist(depvars, 'WRFght'): WRFght_compute = True
220        if gen.searchInlist(depvars, 'WRFdens'): WRFdens_compute = True
221        if gen.searchInlist(depvars, 'WRFpos'): WRFpos_compute = True
222        if gen.searchInlist(depvars, 'WRFtime'): WRFtime_compute = True
223        if gen.searchInlist(depvars, 'WRFz'): WRFz_compute = True
224        if gen.searchInlist(depvars, 'WRFdxdy'): WRFdxdy_compute = True
225        if gen.searchInlist(depvars, 'LONLATdxdy'): LONLATdxdy_compute = True
226
227# Dictionary with the new computed variables to be able to add them
228dictcompvars = {}
229if WRFgeop_compute:
230    print '  ' + main + ': Retrieving geopotential value from WRF as PH + PHB'
231    dimv = ncobj.variables['PH'].shape
232    WRFgeop = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
233
234    # Attributes of the variable
235    Vvals = gen.variables_values('WRFgeop')
236    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
237      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
238
239if WRFp_compute:
240    print '  ' + main + ': Retrieving pressure value from WRF as P + PB'
241    dimv = ncobj.variables['P'].shape
242    WRFp = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
243
244    # Attributes of the variable
245    Vvals = gen.variables_values('WRFp')
246    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
247      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
248
249if WRFght_compute:
250    print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...' 
251    WRFght = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
252
253    # Attributes of the variable
254    Vvals = gen.variables_values('WRFght')
255    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
256      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
257
258if WRFrh_compute:
259    print '    ' + main + ": computing relative humidity from WRF as 'Tetens'" +     \
260      ' equation (T,P) ...'
261    p0=100000.
262    p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
263    tk = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
264    qv = ncobj.variables['QVAPOR'][:]
265
266    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
267    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
268
269    WRFrh = qv/data2
270
271    # Attributes of the variable
272    Vvals = gen.variables_values('WRFrh')
273    dictcompvars['WRFrh'] = {'name': Vvals[0], 'standard_name': Vvals[1],            \
274      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
275
276if WRFt_compute:
277    print '    ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
278    p0=100000.
279    p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
280
281    WRFt = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
282
283    # Attributes of the variable
284    Vvals = gen.variables_values('WRFt')
285    dictcompvars['WRFt'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
286      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
287
288if WRFdens_compute:
289    print '    ' + main + ': computing air density from WRF as ((MU + MUB) * ' +     \
290      'DNW)/g ...'
291
292# Just we need in in absolute values: Size of the central grid cell
293##    dxval = ncobj.getncattr('DX')
294##    dyval = ncobj.getncattr('DY')
295##    mapfac = ncobj.variables['MAPFAC_M'][:]
296##    area = dxval*dyval*mapfac
297
298    mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
299    dnw = ncobj.variables['DNW'][:]
300
301    WRFdens = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]),        \
302      dtype=np.float)
303    levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)
304
305    for it in range(mu.shape[0]):
306        for iz in range(dnw.shape[1]):
307            levval.fill(np.abs(dnw[it,iz]))
308            WRFdens[it,iz,:,:] = levval
309            WRFdens[it,iz,:,:] = mu[it,:,:]*WRFdens[it,iz,:,:]/grav
310
311    # Attributes of the variable
312    Vvals = gen.variables_values('WRFdens')
313    dictcompvars['WRFdens'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
314      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
315
316if WRFpos_compute:
317# WRF positions from the lowest-leftest corner of the matrix
318    print '    ' + main + ': computing position from MAPFAC_M as sqrt(DY*j**2 + ' +  \
319      'DX*x**2)*MAPFAC_M ...'
320
321    mapfac = ncobj.variables['MAPFAC_M'][:]
322
323    distx = np.float(ncobj.getncattr('DX'))
324    disty = np.float(ncobj.getncattr('DY'))
325
326    print 'distx:',distx,'disty:',disty
327
328    dx = mapfac.shape[2]
329    dy = mapfac.shape[1]
330    dt = mapfac.shape[0]
331
332    WRFpos = np.zeros((dt, dy, dx), dtype=np.float)
333
334    for i in range(1,dx):
335        WRFpos[0,0,i] = distx*i/mapfac[0,0,i]
336    for j in range(1,dy):
337        i=0
338        WRFpos[0,j,i] = WRFpos[0,j-1,i] + disty/mapfac[0,j,i]
339        for i in range(1,dx):
340#            WRFpos[0,j,i] = np.sqrt((disty*j)**2. + (distx*i)**2.)/mapfac[0,j,i]
341#            WRFpos[0,j,i] = np.sqrt((disty*j)**2. + (distx*i)**2.)
342             WRFpos[0,j,i] = WRFpos[0,j,i-1] + distx/mapfac[0,j,i]
343
344    for it in range(1,dt):
345        WRFpos[it,:,:] = WRFpos[0,:,:]
346
347if WRFtime_compute:
348    print '    ' + main + ': computing time from WRF as CFtime(Times) ...'
349
350    refdate='19491201000000'
351    tunitsval='minutes'
352
353    timeobj = ncobj.variables['Times']
354    timewrfv = timeobj[:]
355
356    yrref=refdate[0:4]
357    monref=refdate[4:6]
358    dayref=refdate[6:8]
359    horref=refdate[8:10]
360    minref=refdate[10:12]
361    secref=refdate[12:14]
362
363    refdateS = yrref + '-' + monref + '-' + dayref + ' ' + horref + ':' + minref +   \
364      ':' + secref
365
366
367    if len(timeobj.shape) == 2:
368        dt = timeobj.shape[0]
369    else:
370        dt = 1
371    WRFtime = np.zeros((dt), dtype=np.float)
372
373    if len(timeobj.shape) == 2:
374        for it in range(dt):
375            wrfdates = gen.datetimeStr_conversion(timewrfv[it,:],'WRFdatetime',      \
376              'matYmdHMS')
377            WRFtime[it] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
378    else:
379        wrfdates = gen.datetimeStr_conversion(timewrfv[:],'WRFdatetime',             \
380          'matYmdHMS')
381        WRFtime[0] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
382
383    tunits = tunitsval + ' since ' + refdateS
384
385    # Attributes of the variable
386    dictcompvars['WRFtime'] = {'name': 'time', 'standard_name': 'time',              \
387      'long_name': 'time', 'units': tunits, 'calendar': 'gregorian'}
388
389if WRFz_compute:
390    print '  ' + main + ': Retrieving z: height above surface value from WRF as ' +  \
391      'unstagger(PH + PHB)/9.8-hgt'
392    dimv = ncobj.variables['PH'].shape
393    WRFzg = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/9.8
394
395    unzgd = (dimv[0], dimv[1]-1, dimv[2], dimv[3])
396    unzg = np.zeros(unzgd, dtype=np.float)
397    unzg = 0.5*(WRFzg[:,0:dimv[1]-1,:,:] + WRFzg[:,1:dimv[1],:,:])
398
399    WRFz = np.zeros(unzgd, dtype=np.float)
400    for iz in range(dimv[1]-1):
401        WRFz[:,iz,:,:] = unzg[:,iz,:,:] - ncobj.variables['HGT'][:]
402
403    # Attributes of the variable
404    Vvals = gen.variables_values('WRFz')
405    dictcompvars['WRFz'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
406      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
407
408if WRFdxdy_compute:
409    print '  ' + main + ': Retrieving dxdy: real distance between grid points ' +    \
410      'from WRF as dx=(XLONG(i+1)-XLONG(i))*DX/MAPFAC_M, dy=(XLAT(j+1)-XLAT(i))*DY/'+\
411      'MAPFAC_M, ds=sqrt(dx**2+dy**2)'
412    dimv = ncobj.variables['XLONG'].shape
413    WRFlon = ncobj.variables['XLONG'][0,:,:]
414    WRFlat = ncobj.variables['XLAT'][0,:,:]
415    WRFmapfac_m = ncobj.variables['MAPFAC_M'][0,:,:]
416    DX = ncobj.DX
417    DY = ncobj.DY
418
419    dimx = dimv[2]
420    dimy = dimv[1]
421
422    WRFdx = np.zeros((dimy,dimx), dtype=np.float)
423    WRFdy = np.zeros((dimy,dimx), dtype=np.float)
424
425    WRFdx[:,0:dimx-1]=(WRFlon[:,1:dimx]-WRFlon[:,0:dimx-1])*DX/WRFmapfac_m[:,0:dimx-1]
426    WRFdy[0:dimy-1,:]=(WRFlat[1:dimy,:]-WRFlat[0:dimy-1,:])*DY/WRFmapfac_m[0:dimy-1,:]
427    WRFds = np.sqrt(WRFdx**2 + WRFdy**2)
428
429    # Attributes of the variable
430    Vvals = gen.variables_values('WRFdx')
431    dictcompvars['WRFdx'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
432      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
433    Vvals = gen.variables_values('WRFdy')
434    dictcompvars['WRFdy'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
435      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
436    Vvals = gen.variables_values('WRFds')
437    dictcompvars['WRFds'] = {'name': Vvals[0], 'standard_name': Vvals[1],            \
438      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
439
440if LONLATdxdy_compute:
441    print '  ' + main + ': Retrieving dxdy: real distance between grid points ' +    \
442      'from a regular lonlat projection as dx=(lon[i+1]-lon[i])*raddeg*Rearth*' +    \
443      'cos(abs(lat[i])); dy=(lat[j+1]-lat[i])*raddeg*Rearth; ds=sqrt(dx**2+dy**2); '+\
444      'raddeg = pi/180; Rearth=6370.0e03'
445    dimv = ncobj.variables['lon'].shape
446    lon = ncobj.variables['lon'][:]
447    lat = ncobj.variables['lat'][:]
448
449    WRFlon, WRFlat = gen.lonlat2D(lon,lat)
450
451    dimx = WRFlon.shape[1]
452    dimy = WRFlon.shape[0]
453
454    WRFdx = np.zeros((dimy,dimx), dtype=np.float)
455    WRFdy = np.zeros((dimy,dimx), dtype=np.float)
456
457    raddeg = np.pi/180.
458
459    Rearth = fdef.module_definitions.earthradii
460
461    WRFdx[:,0:dimx-1]=(WRFlon[:,1:dimx]-WRFlon[:,0:dimx-1])*raddeg*Rearth*           \
462      np.cos(np.abs(WRFlat[:,0:dimx-1]*raddeg))
463    WRFdy[0:dimy-1,:]=(WRFlat[1:dimy,:]-WRFlat[0:dimy-1,:])*raddeg*Rearth
464    WRFds = np.sqrt(WRFdx**2 + WRFdy**2)
465
466    # Attributes of the variable
467    Vvals = gen.variables_values('WRFdx')
468    dictcompvars['WRFdx'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
469      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
470    Vvals = gen.variables_values('WRFdy')
471    dictcompvars['WRFdy'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
472      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
473    Vvals = gen.variables_values('WRFds')
474    dictcompvars['WRFds'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
475      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
476
477### ## #
478# Going for the diagnostics
479### ## #
480print '  ' + main + ' ...'
481varsadd = []
482
483for idiag in range(Ndiags):
484    print '    diagnostic:',diags[idiag]
485    diagn = diags[idiag].split('|')[0]
486    depvars = diags[idiag].split('|')[1].split('@')
487    if not gen.searchInlist(NONcheckdepvars, diagn):
488        if diags[idiag].split('|')[1].find('@') != -1:
489            depvars = diags[idiag].split('|')[1].split('@')
490            if depvars[0] == 'deaccum': diagn='deaccum'
491            if depvars[0] == 'accum': diagn='accum'
492            for depv in depvars:
493                # Checking without extra arguments of a depending variable (':', separated)
494                if depv.find(':') != -1: depv=depv.split(':')[0]
495                if not ncobj.variables.has_key(depv) and not                         \
496                  gen.searchInlist(NONcheckingvars, depv) and                        \
497                  not gen.searchInlist(methods, depv) and not depvars[0] == 'deaccum'\
498                  and not depvars[0] == 'accum' and not depv[0:2] == 'z=':
499                    print errormsg
500                    print '  ' + main + ": file '" + opts.ncfile +                   \
501                      "' does not have variable '" + depv + "' !!"
502                    quit(-1)
503        else:
504            depvars = diags[idiag].split('|')[1]
505            if not ncobj.variables.has_key(depvars) and not                          \
506              gen.searchInlist(NONcheckingvars, depvars) and                         \
507              not gen.searchInlist(methods, depvars):
508                print errormsg
509                print '  ' + main + ": file '" + opts.ncfile +                       \
510                  "' does not have variable '" + depvars + "' !!"
511                quit(-1)
512
513    print "\n    Computing '" + diagn + "' from: ", depvars, '...'
514
515# acraintot: accumulated total precipitation from WRF RAINC, RAINNC, RAINSH
516    if diagn == 'ACRAINTOT':
517           
518        var0 = ncobj.variables[depvars[0]]
519        var1 = ncobj.variables[depvars[1]]
520        var2 = ncobj.variables[depvars[2]]
521
522        diagout = var0[:] + var1[:] + var2[:]
523
524        dnamesvar = var0.dimensions
525        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
526
527        # Removing the nonChecking variable-dimensions from the initial list
528        varsadd = []
529        for nonvd in NONchkvardims:
530            if gen.searchInlist(dvnamesvar,nonvd): dvnamesvar.remove(nonvd)
531            varsadd.append(nonvd)
532
533        ncvar.insert_variable(ncobj, 'pracc', diagout, dnamesvar, dvnamesvar, newnc)
534
535# accum: acumulation of any variable as (Variable, time [as [tunits]
536#   from/since ....], newvarname)
537    elif diagn == 'accum':
538
539        var0 = ncobj.variables[depvars[0]]
540        var1 = ncobj.variables[depvars[1]]
541
542        dnamesvar = var0.dimensions
543        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
544
545        diagout, diagoutd, diagoutvd = diag.compute_accum(var0,dnamesvar,dvnamesvar)
546        # Removing the nonChecking variable-dimensions from the initial list
547        varsadd = []
548        diagoutvd = list(dvnames)
549        for nonvd in NONchkvardims:
550            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
551            varsadd.append(nonvd)
552
553        CFvarn = ncvar.variables_values(depvars[0])[0]
554
555# Removing the flux
556        if depvars[1] == 'XTIME':
557            dtimeunits = var1.getncattr('description')
558            tunits = dtimeunits.split(' ')[0]
559        else:
560            dtimeunits = var1.getncattr('units')
561            tunits = dtimeunits.split(' ')[0]
562
563        dtime = (var1[1] - var1[0])*diag.timeunits_seconds(tunits)
564
565        ncvar.insert_variable(ncobj, CFvarn + 'acc', diagout*dtime, diagoutd, diagoutvd, newnc)
566
567# cllmh with cldfra, pres
568    elif diagn == 'cllmh':
569           
570        var0 = ncobj.variables[depvars[0]]
571        if depvars[1] == 'WRFp':
572            var1 = WRFp
573        else:
574            var01 = ncobj.variables[depvars[1]]
575            if len(size(var1.shape)) < len(size(var0.shape)):
576                var1 = np.brodcast_arrays(var01,var0)[0]
577            else:
578                var1 = var01
579
580        diagout, diagoutd, diagoutvd = diag.Forcompute_cllmh(var0,var1,dnames,dvnames)
581
582        # Removing the nonChecking variable-dimensions from the initial list
583        varsadd = []
584        for nonvd in NONchkvardims:
585            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
586            varsadd.append(nonvd)
587
588        ncvar.insert_variable(ncobj, 'cll', diagout[0,:], diagoutd, diagoutvd, newnc)
589        ncvar.insert_variable(ncobj, 'clm', diagout[1,:], diagoutd, diagoutvd, newnc)
590        ncvar.insert_variable(ncobj, 'clh', diagout[2,:], diagoutd, diagoutvd, newnc)
591
592# clt with cldfra
593    elif diagn == 'clt':
594           
595        var0 = ncobj.variables[depvars]
596        diagout, diagoutd, diagoutvd = diag.Forcompute_clt(var0,dnames,dvnames)
597
598        # Removing the nonChecking variable-dimensions from the initial list
599        varsadd = []
600        for nonvd in NONchkvardims:
601            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
602            varsadd.append(nonvd)
603           
604        ncvar.insert_variable(ncobj, 'clt', diagout, diagoutd, diagoutvd, newnc)
605
606# convini (pr, time)
607    elif diagn == 'convini':
608           
609        var0 = ncobj.variables[depvars[0]][:]
610        var1 = ncobj.variables[depvars[1]][:]
611        otime = ncobj.variables[depvars[1]]
612
613        dnamesvar = ncobj.variables[depvars[0]].dimensions
614        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
615
616        diagout, diagoutd, diagoutvd  = diag.var_convini(var0, var1, dnames, dvnames)
617
618        ncvar.insert_variable(ncobj, 'convini', diagout, diagoutd, diagoutvd, newnc, \
619          fill=gen.fillValueF)
620        # Getting the right units
621        ovar = newnc.variables['convini']
622        if gen.searchInlist(otime.ncattrs(), 'units'): 
623            tunits = otime.getncattr('units')
624            ncvar.set_attribute(ovar, 'units', tunits)
625            newnc.sync()
626
627# deaccum: deacumulation of any variable as (Variable, time [as [tunits]
628#   from/since ....], newvarname)
629    elif diagn == 'deaccum':
630
631        var0 = ncobj.variables[depvars[0]]
632        var1 = ncobj.variables[depvars[1]]
633
634        dnamesvar = var0.dimensions
635        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
636
637        diagout, diagoutd, diagoutvd = diag.compute_deaccum(var0,dnamesvar,dvnamesvar)
638        # Removing the nonChecking variable-dimensions from the initial list
639        varsadd = []
640        diagoutvd = list(dvnames)
641        for nonvd in NONchkvardims:
642            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
643            varsadd.append(nonvd)
644
645# Transforming to a flux
646        if depvars[1] == 'XTIME':
647            dtimeunits = var1.getncattr('description')
648            tunits = dtimeunits.split(' ')[0]
649        else:
650            dtimeunits = var1.getncattr('units')
651            tunits = dtimeunits.split(' ')[0]
652
653        dtime = (var1[1] - var1[0])*diag.timeunits_seconds(tunits)
654        ncvar.insert_variable(ncobj, depvars[2], diagout/dtime, diagoutd, diagoutvd, \
655          newnc)
656
657# fog_K84: Computation of fog and visibility following Kunkel, (1984) as QCLOUD, QICE
658    elif diagn == 'fog_K84':
659
660        var0 = ncobj.variables[depvars[0]]
661        var1 = ncobj.variables[depvars[1]]
662
663        dnamesvar = list(var0.dimensions)
664        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
665
666        diag1, diag2, diagoutd, diagoutvd = diag.Forcompute_fog_K84(var0, var1,      \
667          dnamesvar, dvnamesvar)
668        # Removing the nonChecking variable-dimensions from the initial list
669        varsadd = []
670        diagoutvd = list(dvnames)
671        for nonvd in NONchkvardims:
672            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
673            varsadd.append(nonvd)
674        ncvar.insert_variable(ncobj, 'fog', diag1, diagoutd, diagoutvd, newnc)
675        ncvar.insert_variable(ncobj, 'fogvisblty', diag2, diagoutd, diagoutvd, newnc)
676
677# fog_RUC: Computation of fog and visibility following Kunkel, (1984) as QVAPOR,
678#    WRFt, WRFp or Q2, T2, PSFC
679    elif diagn == 'fog_RUC':
680
681        var0 = ncobj.variables[depvars[0]]
682        print gen.infmsg
683        if depvars[1] == 'WRFt':
684            print '  ' + main + ": computing '" + diagn + "' using 3D variables !!"
685            var1 = WRFt
686            var2 = WRFp
687        else:
688            print '  ' + main + ": computing '" + diagn + "' using 2D variables !!"
689            var1 = ncobj.variables[depvars[1]]
690            var2 = ncobj.variables[depvars[2]]
691
692        dnamesvar = list(var0.dimensions)
693        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
694
695        diag1, diag2, diagoutd, diagoutvd = diag.Forcompute_fog_RUC(var0, var1, var2,\
696          dnamesvar, dvnamesvar)
697        # Removing the nonChecking variable-dimensions from the initial list
698        varsadd = []
699        diagoutvd = list(dvnames)
700        for nonvd in NONchkvardims:
701            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
702            varsadd.append(nonvd)
703        ncvar.insert_variable(ncobj, 'fog', diag1, diagoutd, diagoutvd, newnc)
704        ncvar.insert_variable(ncobj, 'fogvisblty', diag2, diagoutd, diagoutvd, newnc)
705
706# fog_FRAML50: Computation of fog and visibility following Gultepe, I. and
707#   J.A. Milbrandt, 2010 as QVAPOR, WRFt, WRFp or Q2, T2, PSFC
708    elif diagn == 'fog_FRAML50':
709
710        var0 = ncobj.variables[depvars[0]]
711        print gen.infmsg
712        if depvars[1] == 'WRFt':
713            print '  ' + main + ": computing '" + diagn + "' using 3D variables !!"
714            var1 = WRFt
715            var2 = WRFp
716        else:
717            print '  ' + main + ": computing '" + diagn + "' using 2D variables !!"
718            var1 = ncobj.variables[depvars[1]]
719            var2 = ncobj.variables[depvars[2]]
720
721        dnamesvar = list(var0.dimensions)
722        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
723
724        diag1, diag2, diagoutd, diagoutvd = diag.Forcompute_fog_FRAML50(var0, var1,  \
725          var2, dnamesvar, dvnamesvar)
726        # Removing the nonChecking variable-dimensions from the initial list
727        varsadd = []
728        diagoutvd = list(dvnames)
729        for nonvd in NONchkvardims:
730            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
731            varsadd.append(nonvd)
732        ncvar.insert_variable(ncobj, 'fog', diag1, diagoutd, diagoutvd, newnc)
733        ncvar.insert_variable(ncobj, 'fogvisblty', diag2, diagoutd, diagoutvd, newnc)
734
735# LMDZrh (pres, t, r)
736    elif diagn == 'LMDZrh':
737           
738        var0 = ncobj.variables[depvars[0]][:]
739        var1 = ncobj.variables[depvars[1]][:]
740        var2 = ncobj.variables[depvars[2]][:]
741
742        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnames,dvnames)
743        ncvar.insert_variable(ncobj, 'hur', diagout, diagoutd, diagoutvd, newnc)
744
745# LMDZrhs (psol, t2m, q2m)
746    elif diagn == 'LMDZrhs':
747           
748        var0 = ncobj.variables[depvars[0]][:]
749        var1 = ncobj.variables[depvars[1]][:]
750        var2 = ncobj.variables[depvars[2]][:]
751
752        dnamesvar = ncobj.variables[depvars[0]].dimensions
753        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
754
755        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
756
757        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
758
759# range_faces: LON, LAT, HGT, WRFdxdy, 'face:['WE'/'SN']:[dsfilt]:[dsnewrange]:[hvalleyrange]'
760    elif diagn == 'range_faces':
761           
762        var0 = ncobj.variables[depvars[0]][:]
763        var1 = ncobj.variables[depvars[1]][:]
764        var2 = ncobj.variables[depvars[2]][:]
765        face = depvars[4].split(':')[1]
766        dsfilt = np.float(depvars[4].split(':')[2])
767        dsnewrange = np.float(depvars[4].split(':')[3])
768        hvalleyrange = np.float(depvars[4].split(':')[4])
769
770        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
771        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
772        lon, lat = gen.lonlat2D(var0, var1)
773        if len(var2.shape) == 3:
774            print warnmsg
775            print '  ' + diagn + ": shapping to 2D variable '" + depvars[2] + "' !!"
776            hgt = var2[0,:,:]
777            dnamesvar.pop(0)
778            dvnamesvar.pop(0)           
779        else:
780            hgt = var2[:]
781
782        orogmax, ptorogmax, dhgt, peaks, valleys, origfaces, diagout, diagoutd,      \
783          diagoutvd, rng, rngorogmax, ptrngorogmax= diag.Forcompute_range_faces(lon, \
784           lat, hgt, WRFds, face, dsfilt, dsnewrange, hvalleyrange, dnamesvar,       \
785           dvnamesvar)
786
787        # Removing the nonChecking variable-dimensions from the initial list
788        varsadd = []
789        diagoutvd = list(dvnames)
790        for nonvd in NONchkvardims:
791            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
792            varsadd.append(nonvd)
793
794        ncvar.insert_variable(ncobj, 'dx', WRFdx, diagoutd, diagoutvd, newnc)
795        ncvar.insert_variable(ncobj, 'dy', WRFdy, diagoutd, diagoutvd, newnc)
796        ncvar.insert_variable(ncobj, 'ds', WRFds, diagoutd, diagoutvd, newnc)
797
798        # adding variables to output file
799        if face == 'WE': axis = 'lon'
800        elif face == 'SN': axis = 'lat'
801
802        ncvar.insert_variable(ncobj, 'range', rng, diagoutd, diagoutvd, newnc,       \
803          fill=gen.fillValueI)
804        ovar = newnc.variables['range']
805        ncvar.set_attribute(ovar, 'deriv', axis)
806
807        ncvar.insert_variable(ncobj, 'orogmax', rngorogmax, diagoutd, diagoutvd,     \
808          newnc, fill=gen.fillValueF)
809        newnc.renameVariable('orogmax', 'rangeorogmax')
810        ovar = newnc.variables['rangeorogmax']
811        ncvar.set_attribute(ovar, 'deriv', axis)
812        stdn = ovar.standard_name
813        ncvar.set_attribute(ovar, 'standard_name', 'range_' + stdn)
814        Ln = ovar.long_name
815        ncvar.set_attribute(ovar, 'long_name', 'range ' + stdn)
816
817        ncvar.insert_variable(ncobj, 'ptorogmax', ptrngorogmax, diagoutd, diagoutvd, \
818          newnc)
819        newnc.renameVariable('ptorogmax', 'rangeptorogmax')
820        ovar = newnc.variables['rangeptorogmax']
821        ncvar.set_attribute(ovar, 'deriv', axis)
822        stdn = ovar.standard_name
823        ncvar.set_attribute(ovar, 'standard_name', 'range_' + stdn)
824        Ln = ovar.long_name
825        ncvar.set_attribute(ovar, 'long_name', 'range ' + stdn)
826
827        ncvar.insert_variable(ncobj, 'orogmax', orogmax, diagoutd, diagoutvd,        \
828          newnc)
829        ovar = newnc.variables['orogmax']
830        ncvar.set_attribute(ovar, 'deriv', axis)
831
832        ncvar.insert_variable(ncobj, 'ptorogmax', ptorogmax, diagoutd, diagoutvd,    \
833          newnc)
834        ovar = newnc.variables['ptorogmax']
835        ncvar.set_attribute(ovar, 'deriv', axis)
836
837        ncvar.insert_variable(ncobj, 'orogderiv', dhgt, diagoutd, diagoutvd, newnc)
838        ovar = newnc.variables['orogderiv']
839        ncvar.set_attribute(ovar, 'deriv', axis)
840
841        ncvar.insert_variable(ncobj, 'peak', peaks, diagoutd, diagoutvd, newnc)
842        ncvar.insert_variable(ncobj, 'valley', valleys, diagoutd, diagoutvd, newnc)
843
844        ncvar.insert_variable(ncobj, 'rangefaces', diagout, diagoutd, diagoutvd,    \
845          newnc)
846        newnc.renameVariable('rangefaces', 'rangefacesfilt')
847        ovar = newnc.variables['rangefacesfilt']
848        ncvar.set_attribute(ovar, 'face', face)
849        ncvar.set_attributek(ovar, 'dist_filter', dsfilt, 'F')
850
851        ncvar.insert_variable(ncobj, 'rangefaces', origfaces, diagoutd, diagoutvd,  \
852          newnc, fill=gen.fillValueI)
853        ovar = newnc.variables['rangefaces']
854        ncvar.set_attribute(ovar, 'face', face)
855        ncvar.set_attributek(ovar, 'dist_newrange', dsnewrange, 'F')
856        ncvar.set_attributek(ovar, 'h_valley_newrange', hvalleyrange, 'F')
857
858# mrso: total soil moisture SMOIS, DZS
859    elif diagn == 'WRFmrso':
860           
861        var0 = ncobj.variables[depvars[0]][:]
862        var10 = ncobj.variables[depvars[1]][:]
863        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
864        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
865
866        var1 = var0.copy()*0.
867        var2 = var0.copy()*0.+1.
868        # Must be a better way....
869        for j in range(var0.shape[2]):
870          for i in range(var0.shape[3]):
871              var1[:,:,j,i] = var10
872
873        diagout, diagoutd, diagoutvd = diag.Forcompute_zint(var0, var1, var2,        \
874          dnamesvar, dvnamesvar)
875
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)
882        ncvar.insert_variable(ncobj, 'mrso', diagout, diagoutd, diagoutvd, newnc)
883
884# mrsos: First layer soil moisture SMOIS, DZS
885    elif diagn == 'WRFmrsos':
886           
887        var0 = ncobj.variables[depvars[0]][:]
888        var1 = ncobj.variables[depvars[1]][:]
889        diagoutd = list(ncobj.variables[depvars[0]].dimensions)
890        diagoutvd = ncvar.var_dim_dimv(diagoutd,dnames,dvnames)
891
892        diagoutd.pop(1)
893        diagoutvd.pop(1)
894
895        diagout= np.zeros((var0.shape[0],var0.shape[2],var0.shape[3]), dtype=np.float)
896
897        # Must be a better way....
898        for j in range(var0.shape[2]):
899          for i in range(var0.shape[3]):
900              diagout[:,j,i] = var0[:,0,j,i]*var1[:,0]
901
902        # Removing the nonChecking variable-dimensions from the initial list
903        varsadd = []
904        diagoutvd = list(dvnames)
905        for nonvd in NONchkvardims:
906            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
907            varsadd.append(nonvd)
908        ncvar.insert_variable(ncobj, 'mrsos', diagout, diagoutd, diagoutvd, newnc)
909
910# mslp: mean sea level pressure (pres, psfc, terrain, temp, qv)
911    elif diagn == 'mslp' or diagn == 'WRFmslp':
912           
913        var1 = ncobj.variables[depvars[1]][:]
914        var2 = ncobj.variables[depvars[2]][:]
915        var4 = ncobj.variables[depvars[4]][:]
916
917        if diagn == 'WRFmslp':
918            var0 = WRFp
919            var3 = WRFt
920            dnamesvar = ncobj.variables['P'].dimensions
921        else:
922            var0 = ncobj.variables[depvars[0]][:]
923            var3 = ncobj.variables[depvars[3]][:]
924            dnamesvar = ncobj.variables[depvars[0]].dimensions
925
926        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
927
928        diagout, diagoutd, diagoutvd = diag.compute_mslp(var0, var1, var2, var3, var4,    \
929          dnamesvar, dvnamesvar)
930
931        # Removing the nonChecking variable-dimensions from the initial list
932        varsadd = []
933        diagoutvd = list(dvnames)
934        for nonvd in NONchkvardims:
935            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
936            varsadd.append(nonvd)
937        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
938
939# OMEGAw (omega, p, t) from NCL formulation (https://www.ncl.ucar.edu/Document/Functions/Contributed/omega_to_w.shtml)
940    elif diagn == 'OMEGAw':
941           
942        var0 = ncobj.variables[depvars[0]][:]
943        var1 = ncobj.variables[depvars[1]][:]
944        var2 = ncobj.variables[depvars[2]][:]
945
946        dnamesvar = ncobj.variables[depvars[0]].dimensions
947        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
948
949        diagout, diagoutd, diagoutvd = diag.compute_OMEGAw(var0,var1,var2,dnamesvar,dvnamesvar)
950
951        ncvar.insert_variable(ncobj, 'wa', diagout, diagoutd, diagoutvd, newnc)
952
953# raintot: instantaneous total precipitation from WRF as (RAINC + RAINC + RAINSH) / dTime
954    elif diagn == 'RAINTOT':
955
956        var0 = ncobj.variables[depvars[0]]
957        var1 = ncobj.variables[depvars[1]]
958        var2 = ncobj.variables[depvars[2]]
959
960        if depvars[3] != 'WRFtime':
961            var3 = ncobj.variables[depvars[3]]
962        else:
963            var3 = np.arange(var0.shape[0], dtype=int)
964
965        var = var0[:] + var1[:] + var2[:]
966
967        dnamesvar = var0.dimensions
968        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
969
970        diagout, diagoutd, diagoutvd = diag.compute_deaccum(var,dnamesvar,dvnamesvar)
971
972# Transforming to a flux
973        if var3.shape[0] > 1:
974            if depvars[3] != 'WRFtime':
975                dtimeunits = var3.getncattr('units')
976                tunits = dtimeunits.split(' ')[0]
977   
978                dtime = (var3[1] - var3[0])*diag.timeunits_seconds(tunits)
979            else:
980                var3 = ncobj.variables['Times']
981                time1 = var3[0,:]
982                time2 = var3[1,:]
983                tmf1 = ''
984                tmf2 = ''
985                for ic in range(len(time1)):
986                    tmf1 = tmf1 + time1[ic]
987                    tmf2 = tmf2 + time2[ic]
988                dtdate1 = dtime.datetime.strptime(tmf1,"%Y-%m-%d_%H:%M:%S")
989                dtdate2 = dtime.datetime.strptime(tmf2,"%Y-%m-%d_%H:%M:%S")
990                diffdate12 = dtdate2 - dtdate1
991                dtime = diffdate12.total_seconds()
992                print 'dtime:',dtime
993        else:
994            print warnmsg
995            print '  ' + main + ": only 1 time-step for '" + diag + "' !!"
996            print '    leaving a zero value!'
997            diagout = var0[:]*0.
998            dtime=1.
999
1000        # Removing the nonChecking variable-dimensions from the initial list
1001        varsadd = []
1002        for nonvd in NONchkvardims:
1003            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
1004            varsadd.append(nonvd)
1005           
1006        ncvar.insert_variable(ncobj, 'pr', diagout/dtime, diagoutd, diagoutvd, newnc)
1007
1008# timemax ([varname], time). When a given variable [varname] got its maximum
1009    elif diagn == 'timemax':
1010           
1011        var0 = ncobj.variables[depvars[0]][:]
1012        var1 = ncobj.variables[depvars[1]][:]
1013
1014        otime = ncobj.variables[depvars[1]]
1015
1016        dnamesvar = ncobj.variables[depvars[0]].dimensions
1017        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1018
1019        diagout, diagoutd, diagoutvd  = diag.var_timemax(var0, var1, dnames,         \
1020          dvnames)
1021
1022        ncvar.insert_variable(ncobj, 'timemax', diagout, diagoutd, diagoutvd, newnc, \
1023          fill=gen.fillValueF)
1024        # Getting the right units
1025        ovar = newnc.variables['timemax']
1026        if gen.searchInlist(otime.ncattrs(), 'units'): 
1027            tunits = otime.getncattr('units')
1028            ncvar.set_attribute(ovar, 'units', tunits)
1029            newnc.sync()
1030        ncvar.set_attribute(ovar, 'variable', depvars[0])
1031
1032# timeoverthres ([varname], time, [value], [CFvarn]). When a given variable [varname]   
1033#   overpass a given [value]. Being [CFvarn] the name of the diagnostics in
1034#   variables_values.dat
1035    elif diagn == 'timeoverthres':
1036           
1037        var0 = ncobj.variables[depvars[0]][:]
1038        var1 = ncobj.variables[depvars[1]][:]
1039        var2 = np.float(depvars[2])
1040        var3 = depvars[3]
1041
1042        otime = ncobj.variables[depvars[1]]
1043
1044        dnamesvar = ncobj.variables[depvars[0]].dimensions
1045        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1046
1047        diagout, diagoutd, diagoutvd  = diag.var_timeoverthres(var0, var1, var2,     \
1048          dnames, dvnames)
1049
1050        ncvar.insert_variable(ncobj, var3, diagout, diagoutd, diagoutvd, newnc, \
1051          fill=gen.fillValueF)
1052        # Getting the right units
1053        ovar = newnc.variables[var3]
1054        if gen.searchInlist(otime.ncattrs(), 'units'): 
1055            tunits = otime.getncattr('units')
1056            ncvar.set_attribute(ovar, 'units', tunits)
1057            newnc.sync()
1058        ncvar.set_attribute(ovar, 'overpassed_threshold', var2)
1059
1060# rhs (psfc, t, q) from TimeSeries files
1061    elif diagn == 'TSrhs':
1062           
1063        p0=100000.
1064        var0 = ncobj.variables[depvars[0]][:]
1065        var1 = (ncobj.variables[depvars[1]][:])*(var0/p0)**(2./7.)
1066        var2 = ncobj.variables[depvars[2]][:]
1067
1068        dnamesvar = ncobj.variables[depvars[0]].dimensions
1069        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1070
1071        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
1072
1073        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
1074
1075# slw: total soil liquid water SH2O, DZS
1076    elif diagn == 'WRFslw':
1077           
1078        var0 = ncobj.variables[depvars[0]][:]
1079        var10 = ncobj.variables[depvars[1]][:]
1080        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1081        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1082
1083        var1 = var0.copy()*0.
1084        var2 = var0.copy()*0.+1.
1085        # Must be a better way....
1086        for j in range(var0.shape[2]):
1087          for i in range(var0.shape[3]):
1088              var1[:,:,j,i] = var10
1089
1090        diagout, diagoutd, diagoutvd = diag.Forcompute_zint(var0, var1, var2,        \
1091          dnamesvar, dvnamesvar)
1092
1093        # Removing the nonChecking variable-dimensions from the initial list
1094        varsadd = []
1095        diagoutvd = list(dvnames)
1096        for nonvd in NONchkvardims:
1097            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1098            varsadd.append(nonvd)
1099        ncvar.insert_variable(ncobj, 'slw', diagout, diagoutd, diagoutvd, newnc)
1100
1101# td (psfc, t, q) from TimeSeries files
1102    elif diagn == 'TStd' or diagn == 'td':
1103           
1104        var0 = ncobj.variables[depvars[0]][:]
1105        var1 = ncobj.variables[depvars[1]][:] - 273.15
1106        var2 = ncobj.variables[depvars[2]][:]
1107
1108        dnamesvar = ncobj.variables[depvars[0]].dimensions
1109        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1110
1111        diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
1112
1113        ncvar.insert_variable(ncobj, 'tdas', diagout, diagoutd, diagoutvd, newnc)
1114
1115# td (psfc, t, q) from TimeSeries files
1116    elif diagn == 'TStdC' or diagn == 'tdC':
1117           
1118        var0 = ncobj.variables[depvars[0]][:]
1119# Temperature is already in degrees Celsius
1120        var1 = ncobj.variables[depvars[1]][:]
1121        var2 = ncobj.variables[depvars[2]][:]
1122
1123        dnamesvar = ncobj.variables[depvars[0]].dimensions
1124        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1125
1126        diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
1127
1128        # Removing the nonChecking variable-dimensions from the initial list
1129        varsadd = []
1130        for nonvd in NONchkvardims:
1131            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1132            varsadd.append(nonvd)
1133
1134        ncvar.insert_variable(ncobj, 'tdas', diagout, diagoutd, diagoutvd, newnc)
1135
1136# wds (u, v)
1137    elif diagn == 'TSwds' or diagn == 'wds' :
1138 
1139        var0 = ncobj.variables[depvars[0]][:]
1140        var1 = ncobj.variables[depvars[1]][:]
1141
1142        dnamesvar = ncobj.variables[depvars[0]].dimensions
1143        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1144
1145        diagout, diagoutd, diagoutvd = diag.compute_wds(var0,var1,dnamesvar,dvnamesvar)
1146
1147        # Removing the nonChecking variable-dimensions from the initial list
1148        varsadd = []
1149        for nonvd in NONchkvardims:
1150            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1151            varsadd.append(nonvd)
1152
1153        ncvar.insert_variable(ncobj, 'wds', diagout, diagoutd, diagoutvd, newnc)
1154
1155# wss (u, v)
1156    elif diagn == 'TSwss' or diagn == 'wss':
1157           
1158        var0 = ncobj.variables[depvars[0]][:]
1159        var1 = ncobj.variables[depvars[1]][:]
1160
1161        dnamesvar = ncobj.variables[depvars[0]].dimensions
1162        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1163
1164        diagout, diagoutd, diagoutvd = diag.compute_wss(var0,var1,dnamesvar,dvnamesvar)
1165
1166        # Removing the nonChecking variable-dimensions from the initial list
1167        varsadd = []
1168        for nonvd in NONchkvardims:
1169            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1170            varsadd.append(nonvd)
1171
1172        ncvar.insert_variable(ncobj, 'wss', diagout, diagoutd, diagoutvd, newnc)
1173
1174# turbulence (var)
1175    elif diagn == 'turbulence':
1176
1177        var0 = ncobj.variables[depvars][:]
1178
1179        dnamesvar = list(ncobj.variables[depvars].dimensions)
1180        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1181
1182        diagout, diagoutd, diagoutvd = diag.compute_turbulence(var0,dnamesvar,dvnamesvar)
1183        valsvar = gen.variables_values(depvars)
1184
1185        newvarn = depvars + 'turb'
1186        ncvar.insert_variable(ncobj, newvarn, diagout, diagoutd, 
1187          diagoutvd, newnc)
1188        varobj = newnc.variables[newvarn]
1189        attrv = varobj.long_name
1190        attr = varobj.delncattr('long_name')
1191        newattr = ncvar.set_attribute(varobj, 'long_name', attrv +                   \
1192          " Taylor decomposition turbulence term")
1193
1194# ua va from ws wd (deg)
1195    elif diagn == 'uavaFROMwswd':
1196           
1197        var0 = ncobj.variables[depvars[0]][:]
1198        var1 = ncobj.variables[depvars[1]][:]
1199
1200        ua = var0*np.cos(var1*np.pi/180.)
1201        va = var0*np.sin(var1*np.pi/180.)
1202
1203        dnamesvar = ncobj.variables[depvars[0]].dimensions
1204        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1205
1206        ncvar.insert_variable(ncobj, 'ua', ua, dnamesvar, dvnamesvar, newnc)
1207        ncvar.insert_variable(ncobj, 'va', va, dnamesvar, dvnamesvar, newnc)
1208
1209# ua va from obs ws wd (deg)
1210    elif diagn == 'uavaFROMobswswd':
1211           
1212        var0 = ncobj.variables[depvars[0]][:]
1213        var1 = ncobj.variables[depvars[1]][:]
1214
1215        ua = var0*np.cos((var1+180.)*np.pi/180.)
1216        va = var0*np.sin((var1+180.)*np.pi/180.)
1217
1218        dnamesvar = ncobj.variables[depvars[0]].dimensions
1219        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1220
1221        ncvar.insert_variable(ncobj, 'ua', ua, dnamesvar, dvnamesvar, newnc)
1222        ncvar.insert_variable(ncobj, 'va', va, dnamesvar, dvnamesvar, newnc)
1223
1224# WRFbils fom WRF as HFX + LH
1225    elif diagn == 'WRFbils':
1226           
1227        var0 = ncobj.variables[depvars[0]][:]
1228        var1 = ncobj.variables[depvars[1]][:]
1229
1230        diagout = var0 + var1
1231        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1232        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1233
1234        ncvar.insert_variable(ncobj, 'bils', diagout, dnamesvar, dvnamesvar, newnc)
1235
1236# WRFcape_afwa CAPE, CIN, ZLFC, PLFC, LI following WRF 'phys/module_diaf_afwa.F'
1237#   methodology as WRFt, WRFrh, WRFp, WRFgeop, HGT
1238    elif diagn == 'WRFcape_afwa':
1239        var0 = WRFt
1240        var1 = WRFrh
1241        var2 = WRFp
1242        dz = WRFgeop.shape[1]
1243        # de-staggering
1244        var3 = 0.5*(WRFgeop[:,0:dz-1,:,:]+WRFgeop[:,1:dz,:,:])/9.8
1245        var4 = ncobj.variables[depvars[4]][0,:,:]
1246
1247        dnamesvar = list(ncobj.variables['T'].dimensions)
1248        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1249
1250        diagout = np.zeros(var0.shape, dtype=np.float)
1251        diagout1, diagout2, diagout3, diagout4, diagout5, diagoutd, diagoutvd =      \
1252          diag.Forcompute_cape_afwa(var0, var1, var2, var3, var4, 3, dnamesvar,      \
1253          dvnamesvar)
1254
1255        # Removing the nonChecking variable-dimensions from the initial list
1256        varsadd = []
1257        for nonvd in NONchkvardims:
1258            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1259            varsadd.append(nonvd)
1260
1261        ncvar.insert_variable(ncobj, 'cape', diagout1, diagoutd, diagoutvd, newnc)
1262        ncvar.insert_variable(ncobj, 'cin', diagout2, diagoutd, diagoutvd, newnc)
1263        ncvar.insert_variable(ncobj, 'zlfc', diagout3, diagoutd, diagoutvd, newnc)
1264        ncvar.insert_variable(ncobj, 'plfc', diagout4, diagoutd, diagoutvd, newnc)
1265        ncvar.insert_variable(ncobj, 'li', diagout5, diagoutd, diagoutvd, newnc)
1266
1267# WRFclivi WRF water vapour path WRFdens, QICE, QGRAUPEL, QHAIL
1268    elif diagn == 'WRFclivi':
1269           
1270        var0 = WRFdens
1271        qtot = ncobj.variables[depvars[1]]
1272        qtotv = qtot[:]
1273        Nspecies = len(depvars) - 2
1274        for iv in range(Nspecies):
1275            if ncobj.variables.has_key(depvars[iv+2]):
1276                var1 = ncobj.variables[depvars[iv+2]][:]
1277                qtotv = qtotv + var1
1278
1279        dnamesvar = list(qtot.dimensions)
1280        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1281
1282        diagout, diagoutd, diagoutvd = diag.compute_clivi(var0, qtotv, dnamesvar,dvnamesvar)
1283
1284        # Removing the nonChecking variable-dimensions from the initial list
1285        varsadd = []
1286        diagoutvd = list(dvnames)
1287        for nonvd in NONchkvardims:
1288            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1289            varsadd.append(nonvd)
1290        ncvar.insert_variable(ncobj, 'clivi', diagout, diagoutd, diagoutvd, newnc)
1291
1292# WRFclwvi WRF water cloud-condensed path WRFdens, QCLOUD, QICE, QGRAUPEL, QHAIL
1293    elif diagn == 'WRFclwvi':
1294           
1295        var0 = WRFdens
1296        qtot = ncobj.variables[depvars[1]]
1297        qtotv = ncobj.variables[depvars[1]]
1298        Nspecies = len(depvars) - 2
1299        for iv in range(Nspecies):
1300            if ncobj.variables.has_key(depvars[iv+2]):
1301                var1 = ncobj.variables[depvars[iv+2]]
1302                qtotv = qtotv + var1[:]
1303
1304        dnamesvar = list(qtot.dimensions)
1305        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1306
1307        diagout, diagoutd, diagoutvd = diag.compute_clwvl(var0, qtotv, dnamesvar,dvnamesvar)
1308
1309        # Removing the nonChecking variable-dimensions from the initial list
1310        varsadd = []
1311        diagoutvd = list(dvnames)
1312        for nonvd in NONchkvardims:
1313            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1314            varsadd.append(nonvd)
1315        ncvar.insert_variable(ncobj, 'clwvi', diagout, diagoutd, diagoutvd, newnc)
1316
1317# WRF_denszint WRF vertical integration as WRFdens, sum(Q[water species1], ..., Q[water speciesN]), varn=[varN]
1318    elif diagn == 'WRF_denszint':
1319           
1320        var0 = WRFdens
1321        varn = depvars[1].split('=')[1]
1322        qtot = ncobj.variables[depvars[2]]
1323        qtotv = ncobj.variables[depvars[2]]
1324        Nspecies = len(depvars) - 2
1325        for iv in range(Nspecies):
1326            if ncobj.variables.has_key(depvars[iv+2]):
1327                var1 = ncobj.variables[depvars[iv+2]]
1328                qtotv = qtotv + var1[:]
1329
1330        dnamesvar = list(qtot.dimensions)
1331        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1332
1333        diagout, diagoutd, diagoutvd = diag.compute_clwvl(var0, qtotv, dnamesvar,dvnamesvar)
1334
1335        # Removing the nonChecking variable-dimensions from the initial list
1336        varsadd = []
1337        diagoutvd = list(dvnames)
1338        for nonvd in NONchkvardims:
1339            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1340            varsadd.append(nonvd)
1341        ncvar.insert_variable(ncobj, varn, diagout, diagoutd, diagoutvd, newnc)
1342
1343# WRFgeop geopotential from WRF as PH + PHB
1344    elif diagn == 'WRFgeop':
1345        var0 = ncobj.variables[depvars[0]][:]
1346        var1 = ncobj.variables[depvars[1]][:]
1347
1348        # de-staggering geopotential
1349        diagout0 = var0 + var1
1350        dt = diagout0.shape[0]
1351        dz = diagout0.shape[1]
1352        dy = diagout0.shape[2]
1353        dx = diagout0.shape[3]
1354
1355        diagout = np.zeros((dt,dz-1,dy,dx), dtype=np.float)
1356        diagout = 0.5*(diagout0[:,1:dz,:,:]+diagout0[:,0:dz-1,:,:])
1357
1358        # Removing the nonChecking variable-dimensions from the initial list
1359        varsadd = []
1360        diagoutvd = list(dvnames)
1361        for nonvd in NONchkvardims:
1362            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1363            varsadd.append(nonvd)
1364
1365        ncvar.insert_variable(ncobj, 'zg', diagout, dnames, diagoutvd, newnc)
1366
1367# WRFpotevap_orPM potential evapotranspiration following Penman-Monteith formulation
1368#   implemented in ORCHIDEE (in src_sechiba/enerbil.f90) as: WRFdens, UST, U10, V10, T2, PSFC, QVAPOR
1369    elif diagn == 'WRFpotevap_orPM':
1370        var0 = WRFdens[:,0,:,:]
1371        var1 = ncobj.variables[depvars[1]][:]
1372        var2 = ncobj.variables[depvars[2]][:]
1373        var3 = ncobj.variables[depvars[3]][:]
1374        var4 = ncobj.variables[depvars[4]][:]
1375        var5 = ncobj.variables[depvars[5]][:]
1376        var6 = ncobj.variables[depvars[6]][:,0,:,:]
1377
1378        dnamesvar = list(ncobj.variables[depvars[1]].dimensions)
1379        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1380
1381        diagout = np.zeros(var1.shape, dtype=np.float)
1382        diagout, diagoutd, diagoutvd = diag.Forcompute_potevap_orPM(var0, var1, var2,\
1383          var3, var4, var5, var6, dnamesvar, dvnamesvar)
1384
1385        # Removing the nonChecking variable-dimensions from the initial list
1386        varsadd = []
1387        for nonvd in NONchkvardims:
1388            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1389            varsadd.append(nonvd)
1390
1391        ncvar.insert_variable(ncobj, 'evspsblpot', diagout, diagoutd, diagoutvd, newnc)
1392
1393# WRFmslp_ptarget sea-level pressure following ECMWF method as PSFC, HGT, WRFt, WRFp, ZNU, ZNW
1394    elif diagn == 'WRFpsl_ecmwf':
1395        var0 = ncobj.variables[depvars[0]][:]
1396        var1 = ncobj.variables[depvars[1]][0,:,:]
1397        var2 = WRFt[:,0,:,:]
1398        var4 = WRFp[:,0,:,:]
1399        var5 = ncobj.variables[depvars[4]][0,:]
1400        var6 = ncobj.variables[depvars[5]][0,:]
1401
1402        # This is quite too appriximate!! passing pressure at half-levels to 2nd full
1403        #   level, using eta values at full (ZNW) and half (ZNU) mass levels
1404        var3 = WRFp[:,0,:,:] + (var6[1] - var5[0])*(WRFp[:,1,:,:] - WRFp[:,0,:,:])/  \
1405          (var5[1]-var5[0])
1406
1407        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1408        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1409
1410        diagout = np.zeros(var0.shape, dtype=np.float)
1411        diagout, diagoutd, diagoutvd = diag.Forcompute_psl_ecmwf(var0, var1, var2,   \
1412          var3, var4, dnamesvar, dvnamesvar)
1413
1414        # Removing the nonChecking variable-dimensions from the initial list
1415        varsadd = []
1416        for nonvd in NONchkvardims:
1417            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1418            varsadd.append(nonvd)
1419
1420        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
1421
1422# WRFmslp_ptarget sea-level pressure following ptarget method as WRFp, PSFC, WRFt, HGT, QVAPOR
1423    elif diagn == 'WRFpsl_ptarget':
1424        var0 = WRFp
1425        var1 = ncobj.variables[depvars[1]][:]
1426        var2 = WRFt
1427        var3 = ncobj.variables[depvars[3]][0,:,:]
1428        var4 = ncobj.variables[depvars[4]][:]
1429
1430        dnamesvar = list(ncobj.variables[depvars[4]].dimensions)
1431        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1432
1433        diagout = np.zeros(var0.shape, dtype=np.float)
1434        diagout, diagoutd, diagoutvd = diag.Forcompute_psl_ptarget(var0, var1, var2, \
1435          var3, var4, 700000., dnamesvar, dvnamesvar)
1436
1437        # Removing the nonChecking variable-dimensions from the initial list
1438        varsadd = []
1439        for nonvd in NONchkvardims:
1440            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1441            varsadd.append(nonvd)
1442
1443        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
1444
1445# WRFp pressure from WRF as P + PB
1446    elif diagn == 'WRFp':
1447        var0 = ncobj.variables[depvars[0]][:]
1448        var1 = ncobj.variables[depvars[1]][:]
1449           
1450        diagout = var0 + var1
1451        diagoutd = list(ncobj.variables[depvars[0]].dimensions)
1452        diagoutvd = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1453
1454        # Removing the nonChecking variable-dimensions from the initial list
1455        varsadd = []
1456        for nonvd in NONchkvardims:
1457            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1458            varsadd.append(nonvd)
1459
1460        ncvar.insert_variable(ncobj, 'pres', diagout, diagoutd, diagoutvd, newnc)
1461
1462# WRFpos
1463    elif diagn == 'WRFpos':
1464           
1465        dnamesvar = ncobj.variables['MAPFAC_M'].dimensions
1466        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1467
1468        ncvar.insert_variable(ncobj, 'WRFpos', WRFpos, dnamesvar, dvnamesvar, newnc)
1469
1470# WRFprw WRF water vapour path WRFdens, QVAPOR
1471    elif diagn == 'WRFprw':
1472           
1473        var0 = WRFdens
1474        var1 = ncobj.variables[depvars[1]]
1475
1476        dnamesvar = list(var1.dimensions)
1477        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1478
1479        diagout, diagoutd, diagoutvd = diag.compute_prw(var0, var1, dnamesvar,dvnamesvar)
1480
1481        # Removing the nonChecking variable-dimensions from the initial list
1482        varsadd = []
1483        diagoutvd = list(dvnames)
1484        for nonvd in NONchkvardims:
1485            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1486            varsadd.append(nonvd)
1487        ncvar.insert_variable(ncobj, 'prw', diagout, diagoutd, diagoutvd, newnc)
1488
1489# WRFrh (P, T, QVAPOR)
1490    elif diagn == 'WRFrh':
1491           
1492        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
1493        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1494
1495        ncvar.insert_variable(ncobj, 'hur', WRFrh, dnames, dvnames, newnc)
1496
1497# WRFrhs (PSFC, T2, Q2)
1498    elif diagn == 'WRFrhs':
1499           
1500        var0 = ncobj.variables[depvars[0]][:]
1501        var1 = ncobj.variables[depvars[1]][:]
1502        var2 = ncobj.variables[depvars[2]][:]
1503
1504        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
1505        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1506
1507        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
1508        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
1509
1510# rvors (u10, v10, WRFpos)
1511    elif diagn == 'WRFrvors':
1512           
1513        var0 = ncobj.variables[depvars[0]]
1514        var1 = ncobj.variables[depvars[1]]
1515
1516        diagout = rotational_z(var0, var1, distx)
1517
1518        dnamesvar = ncobj.variables[depvars[0]].dimensions
1519        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1520
1521        ncvar.insert_variable(ncobj, 'rvors', diagout, dnamesvar, dvnamesvar, newnc)
1522
1523# WRFt (T, P, PB)
1524    elif diagn == 'WRFt':
1525        var0 = ncobj.variables[depvars[0]][:]
1526        var1 = ncobj.variables[depvars[1]][:]
1527        var2 = ncobj.variables[depvars[2]][:]
1528
1529        p0=100000.
1530        p=var1 + var2
1531
1532        WRFt = (var0 + 300.)*(p/p0)**(2./7.)
1533
1534        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1535        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1536
1537        # Removing the nonChecking variable-dimensions from the initial list
1538        varsadd = []
1539        diagoutvd = list(dvnames)
1540        for nonvd in NONchkvardims:
1541            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1542            varsadd.append(nonvd)
1543
1544        ncvar.insert_variable(ncobj, 'ta', WRFt, dnames, diagoutvd, newnc)
1545
1546# WRFtda (WRFrh, WRFt)
1547    elif diagn == 'WRFtda':
1548        ARM2 = fdef.module_definitions.arm2
1549        ARM3 = fdef.module_definitions.arm3
1550
1551        gammatarh = np.log(WRFrh) + ARM2*(WRFt-273.15)/((WRFt-273.15)+ARM3)
1552        td = ARM3*gammatarh/(ARM2-gammatarh)
1553
1554        dnamesvar = list(ncobj.variables['T'].dimensions)
1555        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1556
1557        # Removing the nonChecking variable-dimensions from the initial list
1558        varsadd = []
1559        diagoutvd = list(dvnames)
1560        for nonvd in NONchkvardims:
1561            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1562            varsadd.append(nonvd)
1563
1564        ncvar.insert_variable(ncobj, 'tda', td, dnames, diagoutvd, newnc)
1565
1566# WRFtdas (PSFC, T2, Q2)
1567    elif diagn == 'WRFtdas':
1568        ARM2 = fdef.module_definitions.arm2
1569        ARM3 = fdef.module_definitions.arm3
1570
1571        var0 = ncobj.variables[depvars[0]][:]
1572        var1 = ncobj.variables[depvars[1]][:]
1573        var2 = ncobj.variables[depvars[2]][:]
1574
1575        dnamesvar = list(ncobj.variables[depvars[1]].dimensions)
1576        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1577
1578        rhs, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
1579
1580        gammatarhs = np.log(rhs) + ARM2*(var1-273.15)/((var1-273.15)+ARM3)
1581        tdas = ARM3*gammatarhs/(ARM2-gammatarhs) + 273.15
1582
1583        # Removing the nonChecking variable-dimensions from the initial list
1584        varsadd = []
1585        diagoutvd = list(dvnames)
1586        for nonvd in NONchkvardims:
1587            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1588            varsadd.append(nonvd)
1589
1590        ncvar.insert_variable(ncobj, 'tdas', tdas, dnames, diagoutvd, newnc)
1591
1592# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
1593    elif diagn == 'WRFua':
1594        var0 = ncobj.variables[depvars[0]][:]
1595        var1 = ncobj.variables[depvars[1]][:]
1596        var2 = ncobj.variables[depvars[2]][:]
1597        var3 = ncobj.variables[depvars[3]][:]
1598
1599        # un-staggering variables
1600        if len(var0.shape) == 4:
1601            unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1602        elif len(var0.shape) == 3:
1603            # Asuming sunding point (dimt, dimz, dimstgx)
1604            unstgdims = [var0.shape[0], var0.shape[1]]
1605
1606        ua = np.zeros(tuple(unstgdims), dtype=np.float)
1607        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1608        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1609
1610        if len(var0.shape) == 4:
1611            unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1612            unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1613
1614            for iz in range(var0.shape[1]):
1615                ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
1616
1617            dnamesvar = ['Time','bottom_top','south_north','west_east']
1618
1619        elif len(var0.shape) == 3:
1620            unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1])
1621            unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1])
1622            for iz in range(var0.shape[1]):
1623                ua[:,iz] = unstgvar0[:,iz]*var3 - unstgvar1[:,iz]*var2
1624
1625            dnamesvar = ['Time','bottom_top']
1626
1627        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1628
1629        # Removing the nonChecking variable-dimensions from the initial list
1630        varsadd = []
1631        diagoutvd = list(dvnames)
1632        for nonvd in NONchkvardims:
1633            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1634            varsadd.append(nonvd)
1635
1636        ncvar.insert_variable(ncobj, 'ua', ua, dnames, diagoutvd, newnc)
1637
1638# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
1639    elif diagn == 'WRFva':
1640        var0 = ncobj.variables[depvars[0]][:]
1641        var1 = ncobj.variables[depvars[1]][:]
1642        var2 = ncobj.variables[depvars[2]][:]
1643        var3 = ncobj.variables[depvars[3]][:]
1644
1645        # un-staggering variables
1646        if len(var0.shape) == 4:
1647            unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1648        elif len(var0.shape) == 3:
1649            # Asuming sunding point (dimt, dimz, dimstgx)
1650            unstgdims = [var0.shape[0], var0.shape[1]]
1651
1652        va = np.zeros(tuple(unstgdims), dtype=np.float)
1653        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1654        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1655
1656        if len(var0.shape) == 4:
1657            unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1658            unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1659
1660            for iz in range(var0.shape[1]):
1661                va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
1662
1663            dnamesvar = ['Time','bottom_top','south_north','west_east']
1664
1665        elif len(var0.shape) == 3:
1666            unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1])
1667            unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1])
1668            for iz in range(var0.shape[1]):
1669                va[:,iz] = unstgvar0[:,iz]*var2 + unstgvar1[:,iz]*var3
1670
1671            dnamesvar = ['Time','bottom_top']
1672
1673        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1674
1675        # Removing the nonChecking variable-dimensions from the initial list
1676        varsadd = []
1677        diagoutvd = list(dvnames)
1678        for nonvd in NONchkvardims:
1679            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1680            varsadd.append(nonvd)
1681        ncvar.insert_variable(ncobj, 'va', va, dnames, diagoutvd, newnc)
1682
1683
1684# WRFwd (U, V, SINALPHA, COSALPHA) to be rotated !!
1685    elif diagn == 'WRFwd':
1686        var0 = ncobj.variables[depvars[0]][:]
1687        var1 = ncobj.variables[depvars[1]][:]
1688        var2 = ncobj.variables[depvars[2]][:]
1689        var3 = ncobj.variables[depvars[3]][:]
1690
1691        # un-staggering variables
1692        if len(var0.shape) == 4:
1693            unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1694        elif len(var0.shape) == 3:
1695            # Asuming sunding point (dimt, dimz, dimstgx)
1696            unstgdims = [var0.shape[0], var0.shape[1]]
1697
1698        ua = np.zeros(tuple(unstgdims), dtype=np.float)
1699        va = np.zeros(tuple(unstgdims), dtype=np.float)
1700        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1701        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1702
1703        if len(var0.shape) == 4:
1704            unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1705            unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1706
1707            for iz in range(var0.shape[1]):
1708                ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
1709                va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
1710
1711            dnamesvar = ['Time','bottom_top','south_north','west_east']
1712
1713        elif len(var0.shape) == 3:
1714            unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1])
1715            unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1])
1716            for iz in range(var0.shape[1]):
1717                ua[:,iz] = unstgvar0[:,iz]*var3 - unstgvar1[:,iz]*var2
1718                va[:,iz] = unstgvar0[:,iz]*var2 + unstgvar1[:,iz]*var3
1719
1720            dnamesvar = ['Time','bottom_top']
1721
1722        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1723
1724        wd, dnames, dvnames = diag.compute_wd(ua, va, dnamesvar, dvnamesvar)
1725
1726        # Removing the nonChecking variable-dimensions from the initial list
1727        varsadd = []
1728        diagoutvd = list(dvnames)
1729        for nonvd in NONchkvardims:
1730            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1731            varsadd.append(nonvd)
1732
1733        ncvar.insert_variable(ncobj, 'wd', wd, dnames, diagoutvd, newnc)
1734
1735# WRFtime
1736    elif diagn == 'WRFtime':
1737           
1738        diagout = WRFtime
1739
1740        dnamesvar = ['Time']
1741        dvnamesvar = ['Times']
1742
1743        ncvar.insert_variable(ncobj, 'time', diagout, dnamesvar, dvnamesvar, newnc)
1744
1745# ws (U, V)
1746    elif diagn == 'ws':
1747           
1748        # un-staggering variables
1749        if len(var0.shape) == 4:
1750            unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1751        elif len(var0.shape) == 3:
1752            # Asuming sunding point (dimt, dimz, dimstgx)
1753            unstgdims = [var0.shape[0], var0.shape[1]]
1754
1755        ua = np.zeros(tuple(unstgdims), dtype=np.float)
1756        va = np.zeros(tuple(unstgdims), dtype=np.float)
1757        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1758        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1759
1760        if len(var0.shape) == 4:
1761            unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1762            unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1763
1764            for iz in range(var0.shape[1]):
1765                ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
1766                va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
1767
1768            dnamesvar = ['Time','bottom_top','south_north','west_east']
1769
1770        elif len(var0.shape) == 3:
1771            unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1])
1772            unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1])
1773            for iz in range(var0.shape[1]):
1774                ua[:,iz] = unstgvar0[:,iz]*var3 - unstgvar1[:,iz]*var2
1775                va[:,iz] = unstgvar0[:,iz]*var2 + unstgvar1[:,iz]*var3
1776
1777            dnamesvar = ['Time','bottom_top']
1778
1779        diagout = np.sqrt(unstgvar0*unstgvar0 + unstgvar1*unstgvar1)
1780
1781#        dnamesvar = ncobj.variables[depvars[0]].dimensions
1782        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1783
1784        # Removing the nonChecking variable-dimensions from the initial list
1785        varsadd = []
1786        diagoutvd = list(dvnamesvar)
1787        for nonvd in NONchkvardims:
1788            if gen.searchInlist(dvnamesvar,nonvd): diagoutvd.remove(nonvd)
1789            varsadd.append(nonvd)
1790        ncvar.insert_variable(ncobj, 'ws', diagout, dnamesvar, diagoutvd, newnc)
1791
1792# wss (u10, v10)
1793    elif diagn == 'wss':
1794           
1795        var0 = ncobj.variables[depvars[0]][:]
1796        var1 = ncobj.variables[depvars[1]][:]
1797
1798        diagout = np.sqrt(var0*var0 + var1*var1)
1799
1800        dnamesvar = ncobj.variables[depvars[0]].dimensions
1801        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1802
1803        ncvar.insert_variable(ncobj, 'wss', diagout, dnamesvar, dvnamesvar, newnc)
1804
1805# WRFheight height from WRF geopotential as WRFGeop/g
1806    elif diagn == 'WRFheight':
1807           
1808        diagout = WRFgeop/grav
1809
1810        # Removing the nonChecking variable-dimensions from the initial list
1811        varsadd = []
1812        diagoutvd = list(dvnames)
1813        for nonvd in NONchkvardims:
1814            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1815            varsadd.append(nonvd)
1816
1817        ncvar.insert_variable(ncobj, 'zhgt', diagout, dnames, diagoutvd, newnc)
1818
1819# WRFheightrel relative-height from WRF geopotential as WRFgeop(PH + PHB)/g-HGT 'WRFheightrel|PH@PHB@HGT
1820    elif diagn == 'WRFheightrel':
1821        var0 = ncobj.variables[depvars[0]][:]
1822        var1 = ncobj.variables[depvars[1]][:]
1823        var2 = ncobj.variables[depvars[2]][:]
1824
1825        dimz = var0.shape[1]
1826        diagout = np.zeros(tuple(var0.shape), dtype=np.float)
1827        for iz in range(dimz):
1828            diagout[:,iz,:,:] = (var0[:,iz,:,:]+ var1[:,iz,:,:])/grav - var2
1829
1830        # Removing the nonChecking variable-dimensions from the initial list
1831        varsadd = []
1832        diagoutvd = list(dvnames)
1833        for nonvd in NONchkvardims:
1834            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1835            varsadd.append(nonvd)
1836
1837        ncvar.insert_variable(ncobj, 'zhgtrel', diagout, dnames, diagoutvd, newnc)
1838
1839# WRFzmla_gen generic boundary layer hieght computation from WRF theta, QVAPOR, WRFgeop, HGT,
1840    elif diagn == 'WRFzmlagen':
1841        var0 = ncobj.variables[depvars[0]][:]+300.
1842        var1 = ncobj.variables[depvars[1]][:]
1843        dimz = var0.shape[1]
1844        var2 = WRFgeop[:,1:dimz+1,:,:]/9.8
1845        var3 = ncobj.variables[depvars[3]][0,:,:]
1846
1847        diagout, diagoutd, diagoutvd = diag.Forcompute_zmla_gen(var0,var1,var2,var3, \
1848          dnames,dvnames)
1849
1850        # Removing the nonChecking variable-dimensions from the initial list
1851        varsadd = []
1852        for nonvd in NONchkvardims:
1853            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1854            varsadd.append(nonvd)
1855
1856        ncvar.insert_variable(ncobj, 'zmla', diagout, diagoutd, diagoutvd, newnc)
1857
1858# WRFzwind wind extrapolation at a given height using power law computation from WRF
1859#   U, V, WRFz, U10, V10, SINALPHA, COSALPHA, z=[zval]
1860    elif diagn == 'WRFzwind':
1861        var0 = ncobj.variables[depvars[0]][:]
1862        var1 = ncobj.variables[depvars[1]][:]
1863        var2 = WRFz
1864        var3 = ncobj.variables[depvars[3]][:]
1865        var4 = ncobj.variables[depvars[4]][:]
1866        var5 = ncobj.variables[depvars[5]][0,:,:]
1867        var6 = ncobj.variables[depvars[6]][0,:,:]
1868        var7 = np.float(depvars[7].split('=')[1])
1869
1870        # un-staggering 3D winds
1871        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1872        va = np.zeros(tuple(unstgdims), dtype=np.float)
1873        unvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1874        unvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1875        unvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1876        unvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1877
1878        diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwind(unvar0,      \
1879          unvar1, var2, var3, var4, var5, var6, var7, dnames, dvnames)
1880
1881        # Removing the nonChecking variable-dimensions from the initial list
1882        varsadd = []
1883        for nonvd in NONchkvardims:
1884            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1885            varsadd.append(nonvd)
1886
1887        ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
1888        ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)
1889
1890# WRFzwind wind extrapolation at a given hieght using logarithmic law computation
1891#   from WRF U, V, WRFz, U10, V10, SINALPHA, COSALPHA, z=[zval]
1892    elif diagn == 'WRFzwind_log':
1893        var0 = ncobj.variables[depvars[0]][:]
1894        var1 = ncobj.variables[depvars[1]][:]
1895        var2 = WRFz
1896        var3 = ncobj.variables[depvars[3]][:]
1897        var4 = ncobj.variables[depvars[4]][:]
1898        var5 = ncobj.variables[depvars[5]][0,:,:]
1899        var6 = ncobj.variables[depvars[6]][0,:,:]
1900        var7 = np.float(depvars[7].split('=')[1])
1901
1902        # un-staggering 3D winds
1903        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1904        va = np.zeros(tuple(unstgdims), dtype=np.float)
1905        unvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1906        unvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1907        unvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1908        unvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1909
1910        diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwind_log(unvar0,  \
1911          unvar1, var2, var3, var4, var5, var6, var7, dnames, dvnames)
1912
1913        # Removing the nonChecking variable-dimensions from the initial list
1914        varsadd = []
1915        for nonvd in NONchkvardims:
1916            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1917            varsadd.append(nonvd)
1918
1919        ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
1920        ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)
1921
1922# WRFzwindMO wind extrapolation at a given height computation using Monin-Obukhow
1923#   theory from WRF UST, ZNT, RMOL, U10, V10, SINALPHA, COSALPHA, z=[zval]
1924#   NOTE: only useful for [zval] < 80. m
1925    elif diagn == 'WRFzwindMO':
1926        var0 = ncobj.variables[depvars[0]][:]
1927        var1 = ncobj.variables[depvars[1]][:]
1928        var2 = ncobj.variables[depvars[2]][:]
1929        var3 = ncobj.variables[depvars[3]][:]
1930        var4 = ncobj.variables[depvars[4]][:]
1931        var5 = ncobj.variables[depvars[5]][0,:,:]
1932        var6 = ncobj.variables[depvars[6]][0,:,:]
1933        var7 = np.float(depvars[7].split('=')[1])
1934
1935        diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwindMO(var0, var1,\
1936          var2, var3, var4, var5, var6, var7, dnames, dvnames)
1937
1938        # Removing the nonChecking variable-dimensions from the initial list
1939        varsadd = []
1940        for nonvd in NONchkvardims:
1941            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1942            varsadd.append(nonvd)
1943
1944        ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
1945        ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)
1946
1947    else:
1948        print errormsg
1949        print '  ' + main + ": diagnostic '" + diagn + "' not ready!!!"
1950        print '    available diagnostics: ', availdiags
1951        quit(-1)
1952
1953    newnc.sync()
1954    # Adding that additional variables required to compute some diagnostics which
1955    #   where not in the original file
1956    print '  adding additional variables...'
1957    for vadd in varsadd:
1958        if not gen.searchInlist(newnc.variables.keys(),vadd) and                     \
1959          dictcompvars.has_key(vadd):
1960            attrs = dictcompvars[vadd]
1961            vvn = attrs['name']
1962            if not gen.searchInlist(newnc.variables.keys(), vvn):
1963                iidvn = dvnames.index(vadd)
1964                dnn = dnames[iidvn]
1965                if vadd == 'WRFtime':
1966                    dvarvals = WRFtime[:]
1967                newvar = newnc.createVariable(vvn, 'f8', (dnn))
1968                newvar[:] = dvarvals
1969                for attn in attrs.keys():
1970                    if attn != 'name':
1971                        attv = attrs[attn]
1972                        ncvar.set_attribute(newvar, attn, attv)
1973
1974#   end of diagnostics
1975
1976# Global attributes
1977##
1978ncvar.add_global_PyNCplot(newnc, main, None, '2.0')
1979
1980gorigattrs = ncobj.ncattrs()
1981for attr in gorigattrs:
1982    attrv = ncobj.getncattr(attr)
1983    atvar = ncvar.set_attribute(newnc, attr, attrv)
1984
1985ncobj.close()
1986newnc.close()
1987
1988print '\n' + main + ': successfull writting of diagnostics file "' + ofile + '" !!!'
Note: See TracBrowser for help on using the repository browser.