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

Last change on this file since 2656 was 2655, checked in by lfita, 6 years ago

Adding:

  • `front_R04': generalization of front calculation from Rodrigues, 2004
File size: 92.1 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', 'front_R04', 'rhs_tas_tds', 'LMDZrh', 'mslp', 'OMEGAw', 'RAINTOT',      \
93  'range_faces',                                                                     \
94  'rvors', 'td', 'timemax', 'timeoverthres', 'turbulence', 'tws', 'uavaFROMwswd',    \
95  'WRFbnds', 'WRFcape_afwa', 'WRFclivi', 'WRFclwvi', 'WRF_denszint', 'WRFgeop',      \
96  'WRFmrso', 'WRFmrsos', 'WRFpotevap_orPM', 'WRFp', 'WRFpsl_ecmwf',                  \
97  'WRFpsl_ptarget', 'WRFrvors', 'WRFslw', 'ws', 'wds', 'wss', 'WRFheight',           \
98  'WRFheightrel', 'WRFtda', 'WRFtdas', 'WRFtws', 'WRFua', 'WRFva', 'WRFzwind',       \
99  'WRFzwind_log', 'WRFzwindMO', 'zmlagen', 'zmlagenUWsnd']
100
101methods = ['accum', 'deaccum']
102
103# Variables not to check
104NONcheckingvars = ['accum', 'cllmh', 'deaccum', 'face', 'LONLATdxdy',                \
105  'params', 'reglonlatbnds', 'TSrhs', 'TStd', 'TSwds', 'TSwss',                      \
106  'WRFbils',  'WRFbnds',                                                             \
107  'WRFclivi', 'WRFclwvi', 'WRFdens', 'WRFdx', 'WRFdxdy', 'WRFdxdywps', 'WRFdy',      \
108  'WRFgeop', 'WRFp', 'WRFtd',                                                        \
109  'WRFpos', 'WRFprc', 'WRFprls', 'WRFrh', 'LMDZrh', 'LMDZrhs',                       \
110  'WRFrhs', 'WRFrvors',                                                              \
111  'WRFt', 'WRFtime', 'WRFua', 'WRFva', 'WRFwds', 'WRFwss', 'WRFheight', 'WRFz',      \
112  'WRFzg']
113
114# diagnostics not to check their dependeny
115NONcheckdepvars = ['accum', 'deaccum', 'timeoverthres', 'WRF_denszint',              \
116  'WRFzwind_log', 'WRFzwind', 'WRFzwindMO']
117
118NONchkvardims = ['WRFtime']
119
120ofile = 'diagnostics.nc'
121
122dimns = opts.dimns
123varns = opts.varns
124
125# Special method. knowing variable combination
126##
127if opts.dimns == 'variable_combo':
128    print warnmsg
129    print '  ' + main + ': knowing variable combination !!!'
130    combination = variable_combo(opts.varns,opts.ncfile)
131    print '     COMBO: ' + combination
132    quit(-1)
133
134if opts.ncfile is None:
135    print errormsg
136    print '   ' + main + ": No file provided !!"
137    print '     is mandatory to provide a file -f [filename]'
138    quit(-1)
139
140if opts.dimns is None:
141    print errormsg
142    print '   ' + main + ": No description of dimensions are provided !!"
143    print '     is mandatory to provide description of dimensions as -d [dimn]@[vardimname],... '
144    quit(-1)
145
146if opts.varns is None:
147    print errormsg
148    print '   ' + main + ": No variable to diagnose is provided !!"
149    print '     is mandatory to provide a variable to diagnose as -v [diagn]|[varn1]@... '
150    quit(-1)
151
152if not os.path.isfile(opts.ncfile):
153    print errormsg
154    print '   ' + main + ": file '" + opts.ncfile + "' does not exist !!"
155    quit(-1)
156
157ncobj = NetCDFFile(opts.ncfile, 'r')
158
159# Looking for specific variables that might be use in more than one diagnostic
160WRFgeop_compute = False
161WRFp_compute = False
162WRFt_compute = False
163WRFrh_compute = False
164WRFght_compute = False
165WRFdens_compute = False
166WRFpos_compute = False
167WRFtime_compute = False
168WRFz_compute = False
169WRFdx_compute = False
170WRFdx_compute = False
171WRFdxdy_compute = False
172WRFdxdywps_compute = False
173LONLATdxdy_compute = False
174
175# File creation
176newnc = NetCDFFile(ofile,'w')
177
178# dimensions
179dimvalues = dimns.split(',')
180dnames = []
181dvnames = []
182
183for dimval in dimvalues:
184    dn = dimval.split('@')[0]
185    dnv = dimval.split('@')[1]
186    dnames.append(dn)
187    dvnames.append(dnv)
188    # Is there any dimension-variable which should be computed?
189    if dnv == 'WRFgeop':WRFgeop_compute = True
190    if dnv == 'WRFp': WRFp_compute = True
191    if dnv == 'WRFt': WRFt_compute = True
192    if dnv == 'WRFrh': WRFrh_compute = True
193    if dnv == 'WRFght': WRFght_compute = True
194    if dnv == 'WRFdens': WRFdens_compute = True
195    if dnv == 'WRFpos': WRFpos_compute = True
196    if dnv == 'WRFtime': WRFtime_compute = True
197    if dnv == 'WRFz':WRFz_compute = True
198    if dnv == 'WRFdx':WRFdx_compute = True
199    if dnv == 'WRFdy':WRFdy_compute = True
200    if dnv == 'WRFdxdy':WRFdxdy_compute = True
201    if dnv == 'WRFdxdywps':WRFdxdywps_compute = True
202    if dnv == 'LONLATdxdy':LONLATdxdy_compute = True
203
204# diagnostics to compute
205diags = varns.split(',')
206Ndiags = len(diags)
207
208for idiag in range(Ndiags):
209    if diags[idiag].split('|')[1].find('@') == -1:
210        depvars = diags[idiag].split('|')[1]
211        if depvars == 'WRFgeop':WRFgeop_compute = True
212        if depvars == 'WRFp': WRFp_compute = True
213        if depvars == 'WRFt': WRFt_compute = True
214        if depvars == 'WRFrh': WRFrh_compute = True
215        if depvars == 'WRFght': WRFght_compute = True
216        if depvars == 'WRFdens': WRFdens_compute = True
217        if depvars == 'WRFpos': WRFpos_compute = True
218        if depvars == 'WRFtime': WRFtime_compute = True
219        if depvars == 'WRFz': WRFz_compute = True
220    else:
221        depvars = diags[idiag].split('|')[1].split('@')
222        if gen.searchInlist(depvars, 'WRFgeop'): WRFgeop_compute = True
223        if gen.searchInlist(depvars, 'WRFp'): WRFp_compute = True
224        if gen.searchInlist(depvars, 'WRFt'): WRFt_compute = True
225        if gen.searchInlist(depvars, 'WRFrh'): WRFrh_compute = True
226        if gen.searchInlist(depvars, 'WRFght'): WRFght_compute = True
227        if gen.searchInlist(depvars, 'WRFdens'): WRFdens_compute = True
228        if gen.searchInlist(depvars, 'WRFpos'): WRFpos_compute = True
229        if gen.searchInlist(depvars, 'WRFtime'): WRFtime_compute = True
230        if gen.searchInlist(depvars, 'WRFz'): WRFz_compute = True
231        if gen.searchInlist(depvars, 'WRFdx'): WRFdx_compute = True
232        if gen.searchInlist(depvars, 'WRFdy'): WRFdy_compute = True
233        if gen.searchInlist(depvars, 'WRFdxdy'): WRFdxdy_compute = True
234        if gen.searchInlist(depvars, 'WRFdxdywps'): WRFdxdywps_compute = True
235        if gen.searchInlist(depvars, 'LONLATdxdy'): LONLATdxdy_compute = True
236
237# Dictionary with the new computed variables to be able to add them
238dictcompvars = {}
239if WRFgeop_compute:
240    print '  ' + main + ': Retrieving geopotential value from WRF as PH + PHB'
241    dimv = ncobj.variables['PH'].shape
242    WRFgeop = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
243
244    # Attributes of the variable
245    Vvals = gen.variables_values('WRFgeop')
246    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
247      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
248
249if WRFp_compute:
250    print '  ' + main + ': Retrieving pressure value from WRF as P + PB'
251    dimv = ncobj.variables['P'].shape
252    WRFp = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
253
254    # Attributes of the variable
255    Vvals = gen.variables_values('WRFp')
256    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
257      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
258
259if WRFght_compute:
260    print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...' 
261    WRFght = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
262
263    # Attributes of the variable
264    Vvals = gen.variables_values('WRFght')
265    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
266      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
267
268if WRFrh_compute:
269    print '    ' + main + ": computing relative humidity from WRF as 'Tetens'" +     \
270      ' equation (T,P) ...'
271    p0=100000.
272    p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
273    tk = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
274    qv = ncobj.variables['QVAPOR'][:]
275
276    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
277    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
278
279    WRFrh = qv/data2
280
281    # Attributes of the variable
282    Vvals = gen.variables_values('WRFrh')
283    dictcompvars['WRFrh'] = {'name': Vvals[0], 'standard_name': Vvals[1],            \
284      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
285
286if WRFt_compute:
287    print '    ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
288    p0=100000.
289    p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
290
291    WRFt = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
292
293    # Attributes of the variable
294    Vvals = gen.variables_values('WRFt')
295    dictcompvars['WRFt'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
296      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
297
298if WRFdens_compute:
299    print '    ' + main + ': computing air density from WRF as ((MU + MUB) * ' +     \
300      'DNW)/g ...'
301
302# Just we need in in absolute values: Size of the central grid cell
303##    dxval = ncobj.getncattr('DX')
304##    dyval = ncobj.getncattr('DY')
305##    mapfac = ncobj.variables['MAPFAC_M'][:]
306##    area = dxval*dyval*mapfac
307
308    mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
309    dnw = ncobj.variables['DNW'][:]
310
311    WRFdens = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]),        \
312      dtype=np.float)
313    levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)
314
315    for it in range(mu.shape[0]):
316        for iz in range(dnw.shape[1]):
317            levval.fill(np.abs(dnw[it,iz]))
318            WRFdens[it,iz,:,:] = levval
319            WRFdens[it,iz,:,:] = mu[it,:,:]*WRFdens[it,iz,:,:]/grav
320
321    # Attributes of the variable
322    Vvals = gen.variables_values('WRFdens')
323    dictcompvars['WRFdens'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
324      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
325
326if WRFpos_compute:
327# WRF positions from the lowest-leftest corner of the matrix
328    print '    ' + main + ': computing position from MAPFAC_M as sqrt(DY*j**2 + ' +  \
329      'DX*x**2)*MAPFAC_M ...'
330
331    mapfac = ncobj.variables['MAPFAC_M'][:]
332
333    distx = np.float(ncobj.getncattr('DX'))
334    disty = np.float(ncobj.getncattr('DY'))
335
336    print 'distx:',distx,'disty:',disty
337
338    dx = mapfac.shape[2]
339    dy = mapfac.shape[1]
340    dt = mapfac.shape[0]
341
342    WRFpos = np.zeros((dt, dy, dx), dtype=np.float)
343
344    for i in range(1,dx):
345        WRFpos[0,0,i] = distx*i/mapfac[0,0,i]
346    for j in range(1,dy):
347        i=0
348        WRFpos[0,j,i] = WRFpos[0,j-1,i] + disty/mapfac[0,j,i]
349        for i in range(1,dx):
350#            WRFpos[0,j,i] = np.sqrt((disty*j)**2. + (distx*i)**2.)/mapfac[0,j,i]
351#            WRFpos[0,j,i] = np.sqrt((disty*j)**2. + (distx*i)**2.)
352             WRFpos[0,j,i] = WRFpos[0,j,i-1] + distx/mapfac[0,j,i]
353
354    for it in range(1,dt):
355        WRFpos[it,:,:] = WRFpos[0,:,:]
356
357if WRFtime_compute:
358    print '    ' + main + ': computing time from WRF as CFtime(Times) ...'
359
360    refdate='19491201000000'
361    tunitsval='minutes'
362
363    timeobj = ncobj.variables['Times']
364    timewrfv = timeobj[:]
365
366    yrref=refdate[0:4]
367    monref=refdate[4:6]
368    dayref=refdate[6:8]
369    horref=refdate[8:10]
370    minref=refdate[10:12]
371    secref=refdate[12:14]
372
373    refdateS = yrref + '-' + monref + '-' + dayref + ' ' + horref + ':' + minref +   \
374      ':' + secref
375
376
377    if len(timeobj.shape) == 2:
378        dt = timeobj.shape[0]
379    else:
380        dt = 1
381    WRFtime = np.zeros((dt), dtype=np.float)
382
383    if len(timeobj.shape) == 2:
384        for it in range(dt):
385            wrfdates = gen.datetimeStr_conversion(timewrfv[it,:],'WRFdatetime',      \
386              'matYmdHMS')
387            WRFtime[it] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
388    else:
389        wrfdates = gen.datetimeStr_conversion(timewrfv[:],'WRFdatetime',             \
390          'matYmdHMS')
391        WRFtime[0] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
392
393    tunits = tunitsval + ' since ' + refdateS
394
395    # Attributes of the variable
396    dictcompvars['WRFtime'] = {'name': 'time', 'standard_name': 'time',              \
397      'long_name': 'time', 'units': tunits, 'calendar': 'gregorian'}
398
399if WRFz_compute:
400    print '  ' + main + ': Retrieving z: height above surface value from WRF as ' +  \
401      'unstagger(PH + PHB)/9.8-hgt'
402    dimv = ncobj.variables['PH'].shape
403    WRFzg = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/9.8
404
405    unzgd = (dimv[0], dimv[1]-1, dimv[2], dimv[3])
406    unzg = np.zeros(unzgd, dtype=np.float)
407    unzg = 0.5*(WRFzg[:,0:dimv[1]-1,:,:] + WRFzg[:,1:dimv[1],:,:])
408
409    WRFz = np.zeros(unzgd, dtype=np.float)
410    for iz in range(dimv[1]-1):
411        WRFz[:,iz,:,:] = unzg[:,iz,:,:] - ncobj.variables['HGT'][:]
412
413    # Attributes of the variable
414    Vvals = gen.variables_values('WRFz')
415    dictcompvars['WRFz'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
416      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
417
418if WRFdxdy_compute or WRFdx_compute or WRFdy_copute:
419    print '  ' + main + ': Retrieving dxdy: real distance between grid points ' +    \
420      'from WRF as dx=(XLONG(i+1)-XLONG(i))*DX/MAPFAC_M, dy=(XLAT(j+1)-XLAT(i))*DY/'+\
421      'MAPFAC_M, ds=sqrt(dx**2+dy**2)'
422    dimv = ncobj.variables['XLONG'].shape
423    WRFlon = ncobj.variables['XLONG'][0,:,:]
424    WRFlat = ncobj.variables['XLAT'][0,:,:]
425    WRFmapfac_m = ncobj.variables['MAPFAC_M'][0,:,:]
426    DX = ncobj.DX
427    DY = ncobj.DY
428
429    dimx = dimv[2]
430    dimy = dimv[1]
431
432    WRFdx = np.zeros((dimy,dimx), dtype=np.float)
433    WRFdy = np.zeros((dimy,dimx), dtype=np.float)
434
435    WRFdx[:,0:dimx-1]=(WRFlon[:,1:dimx]-WRFlon[:,0:dimx-1])*DX/WRFmapfac_m[:,0:dimx-1]
436    WRFdy[0:dimy-1,:]=(WRFlat[1:dimy,:]-WRFlat[0:dimy-1,:])*DY/WRFmapfac_m[0:dimy-1,:]
437    WRFds = np.sqrt(WRFdx**2 + WRFdy**2)
438
439    # Attributes of the variable
440    Vvals = gen.variables_values('WRFdx')
441    dictcompvars['WRFdx'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
442      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
443    Vvals = gen.variables_values('WRFdy')
444    dictcompvars['WRFdy'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
445      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
446    Vvals = gen.variables_values('WRFds')
447    dictcompvars['WRFds'] = {'name': Vvals[0], 'standard_name': Vvals[1],            \
448      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
449
450if WRFdxdywps_compute:
451    print '  ' + main + ': Retrieving dxdy: real distance between grid points ' +    \
452      'from wpsWRF as dx=(XLONG_M(i+1)-XLONG_M(i))*DX/MAPFAC_M, ' +                  \
453      'dy=(XLAT_M(j+1)-XLAT_M(i))*DY/MAPFAC_M, ds=sqrt(dx**2+dy**2)'
454    dimv = ncobj.variables['XLONG_M'].shape
455    WRFlon = ncobj.variables['XLONG_M'][0,:,:]
456    WRFlat = ncobj.variables['XLAT_M'][0,:,:]
457    WRFmapfac_m = ncobj.variables['MAPFAC_M'][0,:,:]
458    DX = ncobj.DX
459    DY = ncobj.DY
460
461    dimx = dimv[2]
462    dimy = dimv[1]
463
464    WRFdx = np.zeros((dimy,dimx), dtype=np.float)
465    WRFdy = np.zeros((dimy,dimx), dtype=np.float)
466
467    WRFdx[:,0:dimx-1]=(WRFlon[:,1:dimx]-WRFlon[:,0:dimx-1])*DX/WRFmapfac_m[:,0:dimx-1]
468    WRFdy[0:dimy-1,:]=(WRFlat[1:dimy,:]-WRFlat[0:dimy-1,:])*DY/WRFmapfac_m[0:dimy-1,:]
469    WRFds = np.sqrt(WRFdx**2 + WRFdy**2)
470
471    # Attributes of the variable
472    Vvals = gen.variables_values('WRFdx')
473    dictcompvars['WRFdx'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
474      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
475    Vvals = gen.variables_values('WRFdy')
476    dictcompvars['WRFdy'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
477      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
478    Vvals = gen.variables_values('WRFds')
479    dictcompvars['WRFds'] = {'name': Vvals[0], 'standard_name': Vvals[1],            \
480      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
481
482if LONLATdxdy_compute:
483    print '  ' + main + ': Retrieving dxdy: real distance between grid points ' +    \
484      'from a regular lonlat projection as dx=(lon[i+1]-lon[i])*raddeg*Rearth*' +    \
485      'cos(abs(lat[i])); dy=(lat[j+1]-lat[i])*raddeg*Rearth; ds=sqrt(dx**2+dy**2); '+\
486      'raddeg = pi/180; Rearth=6370.0e03'
487    dimv = ncobj.variables['lon'].shape
488    lon = ncobj.variables['lon'][:]
489    lat = ncobj.variables['lat'][:]
490
491    WRFlon, WRFlat = gen.lonlat2D(lon,lat)
492
493    dimx = WRFlon.shape[1]
494    dimy = WRFlon.shape[0]
495
496    WRFdx = np.zeros((dimy,dimx), dtype=np.float)
497    WRFdy = np.zeros((dimy,dimx), dtype=np.float)
498
499    raddeg = np.pi/180.
500
501    Rearth = fdef.module_definitions.earthradii
502
503    WRFdx[:,0:dimx-1]=(WRFlon[:,1:dimx]-WRFlon[:,0:dimx-1])*raddeg*Rearth*           \
504      np.cos(np.abs(WRFlat[:,0:dimx-1]*raddeg))
505    WRFdy[0:dimy-1,:]=(WRFlat[1:dimy,:]-WRFlat[0:dimy-1,:])*raddeg*Rearth
506    WRFds = np.sqrt(WRFdx**2 + WRFdy**2)
507
508    # Attributes of the variable
509    Vvals = gen.variables_values('WRFdx')
510    dictcompvars['WRFdx'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
511      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
512    Vvals = gen.variables_values('WRFdy')
513    dictcompvars['WRFdy'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
514      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
515    Vvals = gen.variables_values('WRFds')
516    dictcompvars['WRFds'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
517      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
518
519### ## #
520# Going for the diagnostics
521### ## #
522print '  ' + main + ' ...'
523varsadd = []
524
525Varns = ncobj.variables.keys()
526Varns.sort()
527
528for idiag in range(Ndiags):
529    print '    diagnostic:',diags[idiag]
530    diagn = diags[idiag].split('|')[0]
531    depvars = diags[idiag].split('|')[1].split('@')
532    if not gen.searchInlist(NONcheckdepvars, diagn):
533        if diags[idiag].split('|')[1].find('@') != -1:
534            depvars = diags[idiag].split('|')[1].split('@')
535            if depvars[0] == 'deaccum': diagn='deaccum'
536            if depvars[0] == 'accum': diagn='accum'
537            for depv in depvars:
538                # Checking without extra arguments of a depending variable (':', separated)
539                if depv.find(':') != -1: depv=depv.split(':')[0]
540                if not ncobj.variables.has_key(depv) and not                         \
541                  gen.searchInlist(NONcheckingvars, depv) and                        \
542                  not gen.searchInlist(methods, depv) and not depvars[0] == 'deaccum'\
543                  and not depvars[0] == 'accum' and not depv[0:2] == 'z=':
544                    print errormsg
545                    print '  ' + main + ": file '" + opts.ncfile +                   \
546                      "' does not have variable '" + depv + "' !!"
547                    print '    available ones:', Varns
548                    quit(-1)
549        else:
550            depvars = diags[idiag].split('|')[1]
551            if not ncobj.variables.has_key(depvars) and not                          \
552              gen.searchInlist(NONcheckingvars, depvars) and                         \
553              not gen.searchInlist(methods, depvars):
554                print errormsg
555                print '  ' + main + ": file '" + opts.ncfile +                       \
556                  "' does not have variable '" + depvars + "' !!"
557                print '    available ones:', Varns
558                quit(-1)
559
560    print "\n    Computing '" + diagn + "' from: ", depvars, '...'
561
562# acraintot: accumulated total precipitation from WRF RAINC, RAINNC, RAINSH
563    if diagn == 'ACRAINTOT':
564           
565        var0 = ncobj.variables[depvars[0]]
566        var1 = ncobj.variables[depvars[1]]
567        var2 = ncobj.variables[depvars[2]]
568
569        diagout = var0[:] + var1[:] + var2[:]
570
571        dnamesvar = var0.dimensions
572        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
573
574        # Removing the nonChecking variable-dimensions from the initial list
575        varsadd = []
576        for nonvd in NONchkvardims:
577            if gen.searchInlist(dvnamesvar,nonvd): dvnamesvar.remove(nonvd)
578            varsadd.append(nonvd)
579
580        ncvar.insert_variable(ncobj, 'pracc', diagout, dnamesvar, dvnamesvar, newnc)
581
582# accum: acumulation of any variable as (Variable, time [as [tunits]
583#   from/since ....], newvarname)
584    elif diagn == 'accum':
585
586        var0 = ncobj.variables[depvars[0]]
587        var1 = ncobj.variables[depvars[1]]
588
589        dnamesvar = var0.dimensions
590        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
591
592        diagout, diagoutd, diagoutvd = diag.compute_accum(var0,dnamesvar,dvnamesvar)
593        # Removing the nonChecking variable-dimensions from the initial list
594        varsadd = []
595        diagoutvd = list(dvnames)
596        for nonvd in NONchkvardims:
597            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
598            varsadd.append(nonvd)
599
600        CFvarn = ncvar.variables_values(depvars[0])[0]
601
602# Removing the flux
603        if depvars[1] == 'XTIME':
604            dtimeunits = var1.getncattr('description')
605            tunits = dtimeunits.split(' ')[0]
606        else:
607            dtimeunits = var1.getncattr('units')
608            tunits = dtimeunits.split(' ')[0]
609
610        dtime = (var1[1] - var1[0])*diag.timeunits_seconds(tunits)
611
612        ncvar.insert_variable(ncobj, CFvarn + 'acc', diagout*dtime, diagoutd, diagoutvd, newnc)
613
614# cllmh with cldfra, pres
615    elif diagn == 'cllmh':
616           
617        var0 = ncobj.variables[depvars[0]]
618        if depvars[1] == 'WRFp':
619            var1 = WRFp
620        else:
621            var01 = ncobj.variables[depvars[1]]
622            if len(size(var1.shape)) < len(size(var0.shape)):
623                var1 = np.brodcast_arrays(var01,var0)[0]
624            else:
625                var1 = var01
626
627        diagout, diagoutd, diagoutvd = diag.Forcompute_cllmh(var0,var1,dnames,dvnames)
628
629        # Removing the nonChecking variable-dimensions from the initial list
630        varsadd = []
631        for nonvd in NONchkvardims:
632            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
633            varsadd.append(nonvd)
634
635        ncvar.insert_variable(ncobj, 'cll', diagout[0,:], diagoutd, diagoutvd, newnc)
636        ncvar.insert_variable(ncobj, 'clm', diagout[1,:], diagoutd, diagoutvd, newnc)
637        ncvar.insert_variable(ncobj, 'clh', diagout[2,:], diagoutd, diagoutvd, newnc)
638
639# clt with cldfra
640    elif diagn == 'clt':
641           
642        var0 = ncobj.variables[depvars]
643        diagout, diagoutd, diagoutvd = diag.Forcompute_clt(var0,dnames,dvnames)
644
645        # Removing the nonChecking variable-dimensions from the initial list
646        varsadd = []
647        for nonvd in NONchkvardims:
648            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
649            varsadd.append(nonvd)
650           
651        ncvar.insert_variable(ncobj, 'clt', diagout, diagoutd, diagoutvd, newnc)
652
653# convini (pr, time)
654    elif diagn == 'convini':
655           
656        var0 = ncobj.variables[depvars[0]][:]
657        var1 = ncobj.variables[depvars[1]][:]
658        otime = ncobj.variables[depvars[1]]
659
660        dnamesvar = ncobj.variables[depvars[0]].dimensions
661        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
662
663        diagout, diagoutd, diagoutvd  = diag.var_convini(var0, var1, dnames, dvnames)
664
665        ncvar.insert_variable(ncobj, 'convini', diagout, diagoutd, diagoutvd, newnc, \
666          fill=gen.fillValueF)
667        # Getting the right units
668        ovar = newnc.variables['convini']
669        if gen.searchInlist(otime.ncattrs(), 'units'): 
670            tunits = otime.getncattr('units')
671            ncvar.set_attribute(ovar, 'units', tunits)
672            newnc.sync()
673
674# deaccum: deacumulation of any variable as (Variable, time [as [tunits]
675#   from/since ....], newvarname)
676    elif diagn == 'deaccum':
677
678        var0 = ncobj.variables[depvars[0]]
679        var1 = ncobj.variables[depvars[1]]
680
681        dnamesvar = var0.dimensions
682        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
683
684        diagout, diagoutd, diagoutvd = diag.compute_deaccum(var0,dnamesvar,dvnamesvar)
685        # Removing the nonChecking variable-dimensions from the initial list
686        varsadd = []
687        diagoutvd = list(dvnames)
688        for nonvd in NONchkvardims:
689            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
690            varsadd.append(nonvd)
691
692# Transforming to a flux
693        if depvars[1] == 'XTIME':
694            dtimeunits = var1.getncattr('description')
695            tunits = dtimeunits.split(' ')[0]
696        else:
697            dtimeunits = var1.getncattr('units')
698            tunits = dtimeunits.split(' ')[0]
699
700        dtime = (var1[1] - var1[0])*diag.timeunits_seconds(tunits)
701        ncvar.insert_variable(ncobj, depvars[2], diagout/dtime, diagoutd, diagoutvd, \
702          newnc)
703
704# fog_K84: Computation of fog and visibility following Kunkel, (1984) as QCLOUD, QICE
705    elif diagn == 'fog_K84':
706
707        var0 = ncobj.variables[depvars[0]]
708        var1 = ncobj.variables[depvars[1]]
709
710        dnamesvar = list(var0.dimensions)
711        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
712
713        diag1, diag2, diagoutd, diagoutvd = diag.Forcompute_fog_K84(var0, var1,      \
714          dnamesvar, dvnamesvar)
715        # Removing the nonChecking variable-dimensions from the initial list
716        varsadd = []
717        diagoutvd = list(dvnames)
718        for nonvd in NONchkvardims:
719            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
720            varsadd.append(nonvd)
721        ncvar.insert_variable(ncobj, 'fog', diag1, diagoutd, diagoutvd, newnc)
722        ncvar.insert_variable(ncobj, 'fogvisblty', diag2, diagoutd, diagoutvd, newnc)
723
724# fog_RUC: Computation of fog and visibility following Kunkel, (1984) as QVAPOR,
725#    WRFt, WRFp or Q2, T2, PSFC
726    elif diagn == 'fog_RUC':
727
728        var0 = ncobj.variables[depvars[0]]
729        print gen.infmsg
730        if depvars[1] == 'WRFt':
731            print '  ' + main + ": computing '" + diagn + "' using 3D variables !!"
732            var1 = WRFt
733            var2 = WRFp
734        else:
735            print '  ' + main + ": computing '" + diagn + "' using 2D variables !!"
736            var1 = ncobj.variables[depvars[1]]
737            var2 = ncobj.variables[depvars[2]]
738
739        dnamesvar = list(var0.dimensions)
740        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
741
742        diag1, diag2, diagoutd, diagoutvd = diag.Forcompute_fog_RUC(var0, var1, var2,\
743          dnamesvar, dvnamesvar)
744        # Removing the nonChecking variable-dimensions from the initial list
745        varsadd = []
746        diagoutvd = list(dvnames)
747        for nonvd in NONchkvardims:
748            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
749            varsadd.append(nonvd)
750        ncvar.insert_variable(ncobj, 'fog', diag1, diagoutd, diagoutvd, newnc)
751        ncvar.insert_variable(ncobj, 'fogvisblty', diag2, diagoutd, diagoutvd, newnc)
752
753# fog_FRAML50: Computation of fog and visibility following Gultepe, I. and
754#   J.A. Milbrandt, 2010 as QVAPOR, WRFt, WRFp or Q2, T2, PSFC
755    elif diagn == 'fog_FRAML50':
756
757        var0 = ncobj.variables[depvars[0]]
758        print gen.infmsg
759        if depvars[1] == 'WRFt':
760            print '  ' + main + ": computing '" + diagn + "' using 3D variables !!"
761            var1 = WRFt
762            var2 = WRFp
763        else:
764            print '  ' + main + ": computing '" + diagn + "' using 2D variables !!"
765            var1 = ncobj.variables[depvars[1]]
766            var2 = ncobj.variables[depvars[2]]
767
768        dnamesvar = list(var0.dimensions)
769        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
770
771        diag1, diag2, diagoutd, diagoutvd = diag.Forcompute_fog_FRAML50(var0, var1,  \
772          var2, dnamesvar, dvnamesvar)
773        # Removing the nonChecking variable-dimensions from the initial list
774        varsadd = []
775        diagoutvd = list(dvnames)
776        for nonvd in NONchkvardims:
777            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
778            varsadd.append(nonvd)
779        ncvar.insert_variable(ncobj, 'fog', diag1, diagoutd, diagoutvd, newnc)
780        ncvar.insert_variable(ncobj, 'fogvisblty', diag2, diagoutd, diagoutvd, newnc)
781
782# front_R04: Computation of the presence of a front as defined by Rodrigues et al.
783#     (2004) (tas, uas, vas, dxs, dys, 'params',[dtas],[dwss])
784#   Rodrigues et al. 2004: Climatologia de frentes frias no litoral de Santa Catarina,
785#     Rev. Bras. Geof. v.22 n.2 Sao Paulo maio/ago. DOI: 10.1590/S0102-261X2004000200004 
786    elif diagn == 'front_R04':
787
788        var0 = ncobj.variables[depvars[0]]
789        var1 = ncobj.variables[depvars[1]]
790        var2 = ncobj.variables[depvars[2]]
791        if depvars[3] == 'WRFdx': var3 = WRFdx
792        else: var3 = ncobj.variables[depvars[3]]
793        if depvars[4] == 'WRFdy': var4 = WRFdy
794        else: var4 = ncobj.variables[depvars[4]]
795        par1 = np.float(depvars[5].split(':')[1])
796        par2 = np.float(depvars[5].split(':')[2])
797
798        dnamesvar = list(var0.dimensions)
799        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
800
801        diag1, diag2, diag3, diag4, diagoutd, diagoutvd = diag.Forcompute_front_R04( \
802          var0, var1, var2, var3, var4, par1, par2, dnamesvar, dvnamesvar)
803
804        # Removing the nonChecking variable-dimensions from the initial list
805        varsadd = []
806        diagoutvd = list(dvnames)
807        for nonvd in NONchkvardims:
808            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
809            varsadd.append(nonvd)
810        ncvar.insert_variable(ncobj, 'front', diag1, diagoutd, diagoutvd, newnc,     \
811          fill=gen.fillValueF)
812        ovar = newnc.variables['front']
813        ovar.setncattr('description', 'front defintion after a generalization of ' + \
814          'Rodrigues et al. 2004, Rev. Bras. Geof. v.22 n.2, ' +                     \
815          'doi: 10.1590/S0102-261X2004000200004')
816        ovar.setncattr('details1', 'N-S wind gradient replaced by gradient of uas, '+\
817          'vas')
818        ovar.setncattr('dtas', par1)
819        ovar.setncattr('dwss', par2)
820        newnc.sync()
821        ncvar.insert_variable(ncobj, 'dt1', diag2, diagoutd, diagoutvd, newnc,       \
822          fill=gen.fillValueF)
823        ovar = newnc.variables['dt1']
824        vu = var0.getncattr('units')
825        ncvar.set_attribute(ovar, 'units', vu)
826        attrv = ovar.getncattr('standard_name')
827        ncvar.set_attribute(ovar, 'standard_name', attrv + depvars[0])
828        attrv = ovar.getncattr('long_name')
829        ncvar.set_attribute(ovar, 'long_name', attrv + ' of ' + depvars[0])
830        ncvar.insert_variable(ncobj, 'gradh', diag3, diagoutd, diagoutvd, newnc,     \
831          fill=gen.fillValueF)
832        ovar = newnc.variables['gradh']
833        vu = var1.getncattr('units')
834        ncvar.set_attribute(ovar, 'units', vu+'m-1')
835        attrv = ovar.getncattr('standard_name')
836        ncvar.set_attribute(ovar, 'standard_name', attrv + 'wss')
837        attrv = ovar.getncattr('long_name')
838        ncvar.set_attribute(ovar, 'long_name', attrv + ' of wss')
839        ncvar.insert_variable(ncobj, 'dt2', diag4, diagoutd, diagoutvd, newnc,       \
840          fill=gen.fillValueF)
841        vu = var0.getncattr('units')
842        ncvar.set_attribute(ovar, 'units', vu)
843        attrv = ovar.getncattr('standard_name')
844        ncvar.set_attribute(ovar, 'standard_name', attrv + depvars[0])
845        attrv = ovar.getncattr('long_name')
846        ncvar.set_attribute(ovar, 'long_name', attrv + ' of ' + depvars[0])
847
848# LMDZrh (pres, t, r)
849    elif diagn == 'LMDZrh':
850           
851        var0 = ncobj.variables[depvars[0]][:]
852        var1 = ncobj.variables[depvars[1]][:]
853        var2 = ncobj.variables[depvars[2]][:]
854
855        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnames,dvnames)
856        ncvar.insert_variable(ncobj, 'hur', diagout, diagoutd, diagoutvd, newnc)
857
858# LMDZrhs (psol, t2m, q2m)
859    elif diagn == 'LMDZrhs':
860           
861        var0 = ncobj.variables[depvars[0]][:]
862        var1 = ncobj.variables[depvars[1]][:]
863        var2 = ncobj.variables[depvars[2]][:]
864
865        dnamesvar = ncobj.variables[depvars[0]].dimensions
866        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
867
868        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
869
870        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
871
872# range_faces: LON, LAT, HGT, WRFdxdy, 'face:['WE'/'SN']:[dsfilt]:[dsnewrange]:[hvalleyrange]'
873    elif diagn == 'range_faces':
874           
875        var0 = ncobj.variables[depvars[0]][:]
876        var1 = ncobj.variables[depvars[1]][:]
877        var2 = ncobj.variables[depvars[2]][:]
878        face = depvars[4].split(':')[1]
879        dsfilt = np.float(depvars[4].split(':')[2])
880        dsnewrange = np.float(depvars[4].split(':')[3])
881        hvalleyrange = np.float(depvars[4].split(':')[4])
882
883        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
884        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
885        lon, lat = gen.lonlat2D(var0, var1)
886        if len(var2.shape) == 3:
887            print warnmsg
888            print '  ' + diagn + ": shapping to 2D variable '" + depvars[2] + "' !!"
889            hgt = var2[0,:,:]
890            dnamesvar.pop(0)
891            dvnamesvar.pop(0)           
892        else:
893            hgt = var2[:]
894
895        orogmax, ptorogmax, dhgt, peaks, valleys, origfaces, diagout, diagoutd,      \
896          diagoutvd, rng, rngorogmax, ptrngorogmax= diag.Forcompute_range_faces(lon, \
897           lat, hgt, WRFdx, WRFdy, WRFds, face, dsfilt, dsnewrange, hvalleyrange,    \
898           dnamesvar, dvnamesvar)
899
900        # Removing the nonChecking variable-dimensions from the initial list
901        varsadd = []
902        diagoutvd = list(dvnames)
903        for nonvd in NONchkvardims:
904            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
905            varsadd.append(nonvd)
906
907        ncvar.insert_variable(ncobj, 'dx', WRFdx, diagoutd, diagoutvd, newnc)
908        ncvar.insert_variable(ncobj, 'dy', WRFdy, diagoutd, diagoutvd, newnc)
909        ncvar.insert_variable(ncobj, 'ds', WRFds, diagoutd, diagoutvd, newnc)
910
911        # adding variables to output file
912        if face == 'WE': axis = 'lon'
913        elif face == 'SN': axis = 'lat'
914
915        ncvar.insert_variable(ncobj, 'range', rng, diagoutd, diagoutvd, newnc,       \
916          fill=gen.fillValueI)
917        ovar = newnc.variables['range']
918        ncvar.set_attribute(ovar, 'deriv', axis)
919
920        ncvar.insert_variable(ncobj, 'orogmax', rngorogmax, diagoutd, diagoutvd,     \
921          newnc, fill=gen.fillValueF)
922        newnc.renameVariable('orogmax', 'rangeorogmax')
923        ovar = newnc.variables['rangeorogmax']
924        ncvar.set_attribute(ovar, 'deriv', axis)
925        stdn = ovar.standard_name
926        ncvar.set_attribute(ovar, 'standard_name', 'range_' + stdn)
927        Ln = ovar.long_name
928        ncvar.set_attribute(ovar, 'long_name', 'range ' + stdn)
929
930        ncvar.insert_variable(ncobj, 'ptorogmax', ptrngorogmax, diagoutd, diagoutvd, \
931          newnc)
932        newnc.renameVariable('ptorogmax', 'rangeptorogmax')
933        ovar = newnc.variables['rangeptorogmax']
934        ncvar.set_attribute(ovar, 'deriv', axis)
935        stdn = ovar.standard_name
936        ncvar.set_attribute(ovar, 'standard_name', 'range_' + stdn)
937        Ln = ovar.long_name
938        ncvar.set_attribute(ovar, 'long_name', 'range ' + stdn)
939
940        ncvar.insert_variable(ncobj, 'orogmax', orogmax, diagoutd, diagoutvd,        \
941          newnc)
942        ovar = newnc.variables['orogmax']
943        ncvar.set_attribute(ovar, 'deriv', axis)
944
945        ncvar.insert_variable(ncobj, 'ptorogmax', ptorogmax, diagoutd, diagoutvd,    \
946          newnc)
947        ovar = newnc.variables['ptorogmax']
948        ncvar.set_attribute(ovar, 'deriv', axis)
949
950        ncvar.insert_variable(ncobj, 'orogderiv', dhgt, diagoutd, diagoutvd, newnc)
951        ovar = newnc.variables['orogderiv']
952        ncvar.set_attribute(ovar, 'deriv', axis)
953
954        ncvar.insert_variable(ncobj, 'peak', peaks, diagoutd, diagoutvd, newnc)
955        ncvar.insert_variable(ncobj, 'valley', valleys, diagoutd, diagoutvd, newnc)
956
957        ncvar.insert_variable(ncobj, 'rangefaces', diagout, diagoutd, diagoutvd,    \
958          newnc)
959        newnc.renameVariable('rangefaces', 'rangefacesfilt')
960        ovar = newnc.variables['rangefacesfilt']
961        ncvar.set_attribute(ovar, 'face', face)
962        ncvar.set_attributek(ovar, 'dist_filter', dsfilt, 'F')
963
964        ncvar.insert_variable(ncobj, 'rangefaces', origfaces, diagoutd, diagoutvd,  \
965          newnc, fill=gen.fillValueI)
966        ovar = newnc.variables['rangefaces']
967        ncvar.set_attribute(ovar, 'face', face)
968        ncvar.set_attributek(ovar, 'dist_newrange', dsnewrange, 'F')
969        ncvar.set_attributek(ovar, 'h_valley_newrange', hvalleyrange, 'F')
970
971# cell_bnds: grid cell bounds from lon, lat from a reglar lon/lat projection  as
972#   intersection of their related parallels and meridians
973    elif diagn == 'reglonlatbnds':
974           
975        var00 = ncobj.variables[depvars[0]][:]
976        var01 = ncobj.variables[depvars[1]][:]
977
978        var0, var1 = gen.lonlat2D(var00,var01)
979
980        dnamesvar = []
981        dnamesvar.append('bnds')
982        if (len(var00.shape) == 3):
983            dnamesvar.append(ncobj.variables[depvars[0]].dimensions[1])
984            dnamesvar.append(ncobj.variables[depvars[0]].dimensions[2])
985        elif (len(var00.shape) == 2):
986            dnamesvar.append(ncobj.variables[depvars[0]].dimensions[0])
987            dnamesvar.append(ncobj.variables[depvars[0]].dimensions[1])
988        elif (len(var00.shape) == 1):
989            dnamesvar.append(ncobj.variables[depvars[0]].dimensions[0])
990            dnamesvar.append(ncobj.variables[depvars[1]].dimensions[0])
991        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
992
993        cellbndsx, cellbndsy, diagoutd, diagoutvd = diag.Forcompute_cellbndsreg(var0,\
994          var1, dnamesvar, dvnamesvar)
995
996        # Removing the nonChecking variable-dimensions from the initial list
997        varsadd = []
998        diagoutvd = list(dvnames)
999        for nonvd in NONchkvardims:
1000            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1001            varsadd.append(nonvd)
1002        # creation of bounds dimension
1003        newdim = newnc.createDimension('bnds', 4)
1004
1005        ncvar.insert_variable(ncobj, 'lon_bnds', cellbndsx, diagoutd, diagoutvd, newnc)
1006        ncvar.insert_variable(ncobj, 'lat_bnds', cellbndsy, diagoutd, diagoutvd, newnc)
1007
1008# tws: Bet Wulb temperature following Stull 2011 (tas, hurs)
1009    elif diagn == 'tws':
1010           
1011        var0 = ncobj.variables[depvars[0]][:]
1012        var1 = ncobj.variables[depvars[1]][:]
1013
1014        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1015        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1016        diagoutd = dnames
1017        diagoutvd = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1018
1019        diagout = diag.var_tws_S11(var0,var1)
1020        ncvar.insert_variable(ncobj, 'tws', diagout, diagoutd, diagoutvd, newnc)
1021
1022# cell_bnds: grid cell bounds from XLONG_U, XLAT_U, XLONG_V, XLAT_V as intersection
1023#   of their related parallels and meridians
1024    elif diagn == 'WRFbnds':
1025           
1026        var0 = ncobj.variables[depvars[0]][0,:,:]
1027        var1 = ncobj.variables[depvars[1]][0,:,:]
1028        var2 = ncobj.variables[depvars[2]][0,:,:]
1029        var3 = ncobj.variables[depvars[3]][0,:,:]
1030
1031        dnamesvar = []
1032        dnamesvar.append('bnds')
1033        dnamesvar.append(ncobj.variables[depvars[0]].dimensions[1])
1034        dnamesvar.append(ncobj.variables[depvars[2]].dimensions[2])
1035        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1036
1037        cellbndsx, cellbndsy, diagoutd, diagoutvd = diag.Forcompute_cellbnds(var0,   \
1038          var1, var2, var3, dnamesvar, dvnamesvar)
1039
1040        # Removing the nonChecking variable-dimensions from the initial list
1041        varsadd = []
1042        diagoutvd = list(dvnames)
1043        for nonvd in NONchkvardims:
1044            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1045            varsadd.append(nonvd)
1046        # creation of bounds dimension
1047        newdim = newnc.createDimension('bnds', 4)
1048
1049        ncvar.insert_variable(ncobj, 'lon_bnds', cellbndsx, diagoutd, diagoutvd, newnc)
1050        newnc.sync()
1051        ncvar.insert_variable(ncobj, 'lat_bnds', cellbndsy, diagoutd, diagoutvd, newnc)
1052        newnc.sync()
1053
1054        if newnc.variables.has_key('XLONG'):
1055            ovar = newnc.variables['XLONG']
1056            ovar.setncattr('bounds', 'lon_bnds lat_bnds')
1057            ovar = newnc.variables['XLAT']
1058            ovar.setncattr('bounds', 'lon_bnds lat_bnds')
1059        elif newnc.variables.has_key('XLONG_M'):
1060            ovar = newnc.variables['XLONG_M']
1061            ovar.setncattr('bounds', 'lon_bnds lat_bnds')
1062            ovar = newnc.variables['XLAT_M']
1063            ovar.setncattr('bounds', 'lon_bnds lat_bnds')
1064        else:
1065            print errormsg
1066            print '  ' + fname + ": error computing diagnostic '" + diagn + "' !!"
1067            print "    neither: 'XLONG', 'XLONG_M' have been found"
1068            quit(-1)
1069
1070# mrso: total soil moisture SMOIS, DZS
1071    elif diagn == 'WRFmrso':
1072           
1073        var0 = ncobj.variables[depvars[0]][:]
1074        var10 = ncobj.variables[depvars[1]][:]
1075        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1076        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1077
1078        var1 = var0.copy()*0.
1079        var2 = var0.copy()*0.+1.
1080        # Must be a better way....
1081        for j in range(var0.shape[2]):
1082          for i in range(var0.shape[3]):
1083              var1[:,:,j,i] = var10
1084
1085        diagout, diagoutd, diagoutvd = diag.Forcompute_zint(var0, var1, var2,        \
1086          dnamesvar, dvnamesvar)
1087
1088        # Removing the nonChecking variable-dimensions from the initial list
1089        varsadd = []
1090        diagoutvd = list(dvnames)
1091        for nonvd in NONchkvardims:
1092            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1093            varsadd.append(nonvd)
1094        ncvar.insert_variable(ncobj, 'mrso', diagout, diagoutd, diagoutvd, newnc)
1095
1096# mrsos: First layer soil moisture SMOIS, DZS
1097    elif diagn == 'WRFmrsos':
1098           
1099        var0 = ncobj.variables[depvars[0]][:]
1100        var1 = ncobj.variables[depvars[1]][:]
1101        diagoutd = list(ncobj.variables[depvars[0]].dimensions)
1102        diagoutvd = ncvar.var_dim_dimv(diagoutd,dnames,dvnames)
1103
1104        diagoutd.pop(1)
1105        diagoutvd.pop(1)
1106
1107        diagout= np.zeros((var0.shape[0],var0.shape[2],var0.shape[3]), dtype=np.float)
1108
1109        # Must be a better way....
1110        for j in range(var0.shape[2]):
1111          for i in range(var0.shape[3]):
1112              diagout[:,j,i] = var0[:,0,j,i]*var1[:,0]
1113
1114        # Removing the nonChecking variable-dimensions from the initial list
1115        varsadd = []
1116        diagoutvd = list(dvnames)
1117        for nonvd in NONchkvardims:
1118            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1119            varsadd.append(nonvd)
1120        ncvar.insert_variable(ncobj, 'mrsos', diagout, diagoutd, diagoutvd, newnc)
1121
1122# mslp: mean sea level pressure (pres, psfc, terrain, temp, qv)
1123    elif diagn == 'mslp' or diagn == 'WRFmslp':
1124           
1125        var1 = ncobj.variables[depvars[1]][:]
1126        var2 = ncobj.variables[depvars[2]][:]
1127        var4 = ncobj.variables[depvars[4]][:]
1128
1129        if diagn == 'WRFmslp':
1130            var0 = WRFp
1131            var3 = WRFt
1132            dnamesvar = ncobj.variables['P'].dimensions
1133        else:
1134            var0 = ncobj.variables[depvars[0]][:]
1135            var3 = ncobj.variables[depvars[3]][:]
1136            dnamesvar = ncobj.variables[depvars[0]].dimensions
1137
1138        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1139
1140        diagout, diagoutd, diagoutvd = diag.compute_mslp(var0, var1, var2, var3, var4,    \
1141          dnamesvar, dvnamesvar)
1142
1143        # Removing the nonChecking variable-dimensions from the initial list
1144        varsadd = []
1145        diagoutvd = list(dvnames)
1146        for nonvd in NONchkvardims:
1147            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1148            varsadd.append(nonvd)
1149        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
1150
1151# WRFtws: Bet Wulb temperature following Stull 2011 (PSFC, T2, Q2)
1152    elif diagn == 'WRFtws':
1153
1154        var0 = ncobj.variables[depvars[0]][:]
1155        var1 = ncobj.variables[depvars[1]][:]
1156        var2 = ncobj.variables[depvars[2]][:]
1157
1158        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
1159        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1160
1161        hurs, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
1162           
1163        diagout = diag.var_tws_S11(var1, hurs)
1164        # Removing the nonChecking variable-dimensions from the initial list
1165        varsadd = []
1166        diagoutvd = list(dvnames)
1167        for nonvd in NONchkvardims:
1168            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1169            varsadd.append(nonvd)
1170
1171        ncvar.insert_variable(ncobj, 'tws', diagout, diagoutd, diagoutvd, newnc)
1172
1173# OMEGAw (omega, p, t) from NCL formulation (https://www.ncl.ucar.edu/Document/Functions/Contributed/omega_to_w.shtml)
1174    elif diagn == 'OMEGAw':
1175           
1176        var0 = ncobj.variables[depvars[0]][:]
1177        var1 = ncobj.variables[depvars[1]][:]
1178        var2 = ncobj.variables[depvars[2]][:]
1179
1180        dnamesvar = ncobj.variables[depvars[0]].dimensions
1181        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1182
1183        diagout, diagoutd, diagoutvd = diag.compute_OMEGAw(var0,var1,var2,dnamesvar,dvnamesvar)
1184
1185        ncvar.insert_variable(ncobj, 'wa', diagout, diagoutd, diagoutvd, newnc)
1186
1187# raintot: instantaneous total precipitation from WRF as (RAINC + RAINC + RAINSH) / dTime
1188    elif diagn == 'RAINTOT':
1189
1190        var0 = ncobj.variables[depvars[0]]
1191        var1 = ncobj.variables[depvars[1]]
1192        var2 = ncobj.variables[depvars[2]]
1193
1194        if depvars[3] != 'WRFtime':
1195            var3 = ncobj.variables[depvars[3]]
1196        else:
1197            var3 = np.arange(var0.shape[0], dtype=int)
1198
1199        var = var0[:] + var1[:] + var2[:]
1200
1201        dnamesvar = var0.dimensions
1202        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1203
1204        diagout, diagoutd, diagoutvd = diag.compute_deaccum(var,dnamesvar,dvnamesvar)
1205
1206# Transforming to a flux
1207        if var3.shape[0] > 1:
1208            if depvars[3] != 'WRFtime':
1209                dtimeunits = var3.getncattr('units')
1210                tunits = dtimeunits.split(' ')[0]
1211   
1212                dtime = (var3[1] - var3[0])*diag.timeunits_seconds(tunits)
1213            else:
1214                var3 = ncobj.variables['Times']
1215                time1 = var3[0,:]
1216                time2 = var3[1,:]
1217                tmf1 = ''
1218                tmf2 = ''
1219                for ic in range(len(time1)):
1220                    tmf1 = tmf1 + time1[ic]
1221                    tmf2 = tmf2 + time2[ic]
1222                dtdate1 = dtime.datetime.strptime(tmf1,"%Y-%m-%d_%H:%M:%S")
1223                dtdate2 = dtime.datetime.strptime(tmf2,"%Y-%m-%d_%H:%M:%S")
1224                diffdate12 = dtdate2 - dtdate1
1225                dtime = diffdate12.total_seconds()
1226                print 'dtime:',dtime
1227        else:
1228            print warnmsg
1229            print '  ' + main + ": only 1 time-step for '" + diag + "' !!"
1230            print '    leaving a zero value!'
1231            diagout = var0[:]*0.
1232            dtime=1.
1233
1234        # Removing the nonChecking variable-dimensions from the initial list
1235        varsadd = []
1236        for nonvd in NONchkvardims:
1237            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
1238            varsadd.append(nonvd)
1239           
1240        ncvar.insert_variable(ncobj, 'pr', diagout/dtime, diagoutd, diagoutvd, newnc)
1241
1242# timemax ([varname], time). When a given variable [varname] got its maximum
1243    elif diagn == 'timemax':
1244           
1245        var0 = ncobj.variables[depvars[0]][:]
1246        var1 = ncobj.variables[depvars[1]][:]
1247
1248        otime = ncobj.variables[depvars[1]]
1249
1250        dnamesvar = ncobj.variables[depvars[0]].dimensions
1251        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1252
1253        diagout, diagoutd, diagoutvd  = diag.var_timemax(var0, var1, dnames,         \
1254          dvnames)
1255
1256        ncvar.insert_variable(ncobj, 'timemax', diagout, diagoutd, diagoutvd, newnc, \
1257          fill=gen.fillValueF)
1258        # Getting the right units
1259        ovar = newnc.variables['timemax']
1260        if gen.searchInlist(otime.ncattrs(), 'units'): 
1261            tunits = otime.getncattr('units')
1262            ncvar.set_attribute(ovar, 'units', tunits)
1263            newnc.sync()
1264        ncvar.set_attribute(ovar, 'variable', depvars[0])
1265
1266# timeoverthres ([varname], time, [value], [CFvarn]). When a given variable [varname]   
1267#   overpass a given [value]. Being [CFvarn] the name of the diagnostics in
1268#   variables_values.dat
1269    elif diagn == 'timeoverthres':
1270           
1271        var0 = ncobj.variables[depvars[0]][:]
1272        var1 = ncobj.variables[depvars[1]][:]
1273        var2 = np.float(depvars[2])
1274        var3 = depvars[3]
1275
1276        otime = ncobj.variables[depvars[1]]
1277
1278        dnamesvar = ncobj.variables[depvars[0]].dimensions
1279        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1280
1281        diagout, diagoutd, diagoutvd  = diag.var_timeoverthres(var0, var1, var2,     \
1282          dnames, dvnames)
1283
1284        ncvar.insert_variable(ncobj, var3, diagout, diagoutd, diagoutvd, newnc, \
1285          fill=gen.fillValueF)
1286        # Getting the right units
1287        ovar = newnc.variables[var3]
1288        if gen.searchInlist(otime.ncattrs(), 'units'): 
1289            tunits = otime.getncattr('units')
1290            ncvar.set_attribute(ovar, 'units', tunits)
1291            newnc.sync()
1292        ncvar.set_attribute(ovar, 'overpassed_threshold', var2)
1293
1294# rhs (psfc, t, q) from TimeSeries files
1295    elif diagn == 'TSrhs':
1296           
1297        p0=100000.
1298        var0 = ncobj.variables[depvars[0]][:]
1299        var1 = (ncobj.variables[depvars[1]][:])*(var0/p0)**(2./7.)
1300        var2 = ncobj.variables[depvars[2]][:]
1301
1302        dnamesvar = ncobj.variables[depvars[0]].dimensions
1303        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1304
1305        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
1306
1307        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
1308
1309# rhs (psfc, t, q) from tas, tds
1310    elif diagn == 'rhs_tas_tds':
1311           
1312        var0 = ncobj.variables[depvars[0]][:]
1313        var1 = ncobj.variables[depvars[1]][:]
1314
1315        dnamesvar = ncobj.variables[depvars[0]].dimensions
1316        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1317
1318        diagout, diagoutd, diagoutvd = diag.var_hur_tas_tds(var0,var1,dnamesvar,     \
1319          dvnamesvar)
1320
1321        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
1322
1323# slw: total soil liquid water SH2O, DZS
1324    elif diagn == 'WRFslw':
1325           
1326        var0 = ncobj.variables[depvars[0]][:]
1327        var10 = ncobj.variables[depvars[1]][:]
1328        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1329        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1330
1331        var1 = var0.copy()*0.
1332        var2 = var0.copy()*0.+1.
1333        # Must be a better way....
1334        for j in range(var0.shape[2]):
1335          for i in range(var0.shape[3]):
1336              var1[:,:,j,i] = var10
1337
1338        diagout, diagoutd, diagoutvd = diag.Forcompute_zint(var0, var1, var2,        \
1339          dnamesvar, dvnamesvar)
1340
1341        # Removing the nonChecking variable-dimensions from the initial list
1342        varsadd = []
1343        diagoutvd = list(dvnames)
1344        for nonvd in NONchkvardims:
1345            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1346            varsadd.append(nonvd)
1347        ncvar.insert_variable(ncobj, 'slw', diagout, diagoutd, diagoutvd, newnc)
1348
1349# td (psfc, t, q) from TimeSeries files
1350    elif diagn == 'TStd' or diagn == 'td':
1351           
1352        var0 = ncobj.variables[depvars[0]][:]
1353        var1 = ncobj.variables[depvars[1]][:] - 273.15
1354        var2 = ncobj.variables[depvars[2]][:]
1355
1356        dnamesvar = ncobj.variables[depvars[0]].dimensions
1357        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1358
1359        diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
1360
1361        ncvar.insert_variable(ncobj, 'tdas', diagout, diagoutd, diagoutvd, newnc)
1362
1363# td (psfc, t, q) from TimeSeries files
1364    elif diagn == 'TStdC' or diagn == 'tdC':
1365           
1366        var0 = ncobj.variables[depvars[0]][:]
1367# Temperature is already in degrees Celsius
1368        var1 = ncobj.variables[depvars[1]][:]
1369        var2 = ncobj.variables[depvars[2]][:]
1370
1371        dnamesvar = ncobj.variables[depvars[0]].dimensions
1372        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1373
1374        diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
1375
1376        # Removing the nonChecking variable-dimensions from the initial list
1377        varsadd = []
1378        for nonvd in NONchkvardims:
1379            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1380            varsadd.append(nonvd)
1381
1382        ncvar.insert_variable(ncobj, 'tdas', diagout, diagoutd, diagoutvd, newnc)
1383
1384# wds (u, v)
1385    elif diagn == 'TSwds' or diagn == 'wds' :
1386 
1387        var0 = ncobj.variables[depvars[0]][:]
1388        var1 = ncobj.variables[depvars[1]][:]
1389
1390        dnamesvar = ncobj.variables[depvars[0]].dimensions
1391        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1392
1393        diagout, diagoutd, diagoutvd = diag.compute_wds(var0,var1,dnamesvar,dvnamesvar)
1394
1395        # Removing the nonChecking variable-dimensions from the initial list
1396        varsadd = []
1397        for nonvd in NONchkvardims:
1398            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1399            varsadd.append(nonvd)
1400
1401        ncvar.insert_variable(ncobj, 'wds', diagout, diagoutd, diagoutvd, newnc)
1402
1403# wss (u, v)
1404    elif diagn == 'TSwss' or diagn == 'wss':
1405           
1406        var0 = ncobj.variables[depvars[0]][:]
1407        var1 = ncobj.variables[depvars[1]][:]
1408
1409        dnamesvar = ncobj.variables[depvars[0]].dimensions
1410        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1411
1412        diagout, diagoutd, diagoutvd = diag.compute_wss(var0,var1,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, 'wss', diagout, diagoutd, diagoutvd, newnc)
1421
1422# turbulence (var)
1423    elif diagn == 'turbulence':
1424
1425        var0 = ncobj.variables[depvars][:]
1426
1427        dnamesvar = list(ncobj.variables[depvars].dimensions)
1428        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1429
1430        diagout, diagoutd, diagoutvd = diag.compute_turbulence(var0,dnamesvar,dvnamesvar)
1431        valsvar = gen.variables_values(depvars)
1432
1433        newvarn = depvars + 'turb'
1434        ncvar.insert_variable(ncobj, newvarn, diagout, diagoutd, 
1435          diagoutvd, newnc)
1436        varobj = newnc.variables[newvarn]
1437        attrv = varobj.long_name
1438        attr = varobj.delncattr('long_name')
1439        newattr = ncvar.set_attribute(varobj, 'long_name', attrv +                   \
1440          " Taylor decomposition turbulence term")
1441
1442# ua va from ws wd (deg)
1443    elif diagn == 'uavaFROMwswd':
1444           
1445        var0 = ncobj.variables[depvars[0]][:]
1446        var1 = ncobj.variables[depvars[1]][:]
1447
1448        ua = var0*np.cos(var1*np.pi/180.)
1449        va = var0*np.sin(var1*np.pi/180.)
1450
1451        dnamesvar = ncobj.variables[depvars[0]].dimensions
1452        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1453
1454        ncvar.insert_variable(ncobj, 'ua', ua, dnamesvar, dvnamesvar, newnc)
1455        ncvar.insert_variable(ncobj, 'va', va, dnamesvar, dvnamesvar, newnc)
1456
1457# ua va from obs ws wd (deg)
1458    elif diagn == 'uavaFROMobswswd':
1459           
1460        var0 = ncobj.variables[depvars[0]][:]
1461        var1 = ncobj.variables[depvars[1]][:]
1462
1463        ua = var0*np.cos((var1+180.)*np.pi/180.)
1464        va = var0*np.sin((var1+180.)*np.pi/180.)
1465
1466        dnamesvar = ncobj.variables[depvars[0]].dimensions
1467        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1468
1469        ncvar.insert_variable(ncobj, 'ua', ua, dnamesvar, dvnamesvar, newnc)
1470        ncvar.insert_variable(ncobj, 'va', va, dnamesvar, dvnamesvar, newnc)
1471
1472# WRFbils fom WRF as HFX + LH
1473    elif diagn == 'WRFbils':
1474           
1475        var0 = ncobj.variables[depvars[0]][:]
1476        var1 = ncobj.variables[depvars[1]][:]
1477
1478        diagout = var0 + var1
1479        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1480        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1481
1482        ncvar.insert_variable(ncobj, 'bils', diagout, dnamesvar, dvnamesvar, newnc)
1483
1484# WRFcape_afwa CAPE, CIN, ZLFC, PLFC, LI following WRF 'phys/module_diaf_afwa.F'
1485#   methodology as WRFt, WRFrh, WRFp, WRFgeop, HGT
1486    elif diagn == 'WRFcape_afwa':
1487        var0 = WRFt
1488        var1 = WRFrh
1489        var2 = WRFp
1490        dz = WRFgeop.shape[1]
1491        # de-staggering
1492        var3 = 0.5*(WRFgeop[:,0:dz-1,:,:]+WRFgeop[:,1:dz,:,:])/9.8
1493        var4 = ncobj.variables[depvars[4]][0,:,:]
1494
1495        dnamesvar = list(ncobj.variables['T'].dimensions)
1496        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1497
1498        diagout = np.zeros(var0.shape, dtype=np.float)
1499        diagout1, diagout2, diagout3, diagout4, diagout5, diagoutd, diagoutvd =      \
1500          diag.Forcompute_cape_afwa(var0, var1, var2, var3, var4, 3, dnamesvar,      \
1501          dvnamesvar)
1502
1503        # Removing the nonChecking variable-dimensions from the initial list
1504        varsadd = []
1505        for nonvd in NONchkvardims:
1506            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1507            varsadd.append(nonvd)
1508
1509        ncvar.insert_variable(ncobj, 'cape', diagout1, diagoutd, diagoutvd, newnc)
1510        ncvar.insert_variable(ncobj, 'cin', diagout2, diagoutd, diagoutvd, newnc)
1511        ncvar.insert_variable(ncobj, 'zlfc', diagout3, diagoutd, diagoutvd, newnc)
1512        ncvar.insert_variable(ncobj, 'plfc', diagout4, diagoutd, diagoutvd, newnc)
1513        ncvar.insert_variable(ncobj, 'li', diagout5, diagoutd, diagoutvd, newnc)
1514
1515# WRFclivi WRF water vapour path WRFdens, QICE, QGRAUPEL, QHAIL
1516    elif diagn == 'WRFclivi':
1517           
1518        var0 = WRFdens
1519        qtot = ncobj.variables[depvars[1]]
1520        qtotv = qtot[:]
1521        Nspecies = len(depvars) - 2
1522        for iv in range(Nspecies):
1523            if ncobj.variables.has_key(depvars[iv+2]):
1524                var1 = ncobj.variables[depvars[iv+2]][:]
1525                qtotv = qtotv + var1
1526
1527        dnamesvar = list(qtot.dimensions)
1528        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1529
1530        diagout, diagoutd, diagoutvd = diag.compute_clivi(var0, qtotv, dnamesvar,dvnamesvar)
1531
1532        # Removing the nonChecking variable-dimensions from the initial list
1533        varsadd = []
1534        diagoutvd = list(dvnames)
1535        for nonvd in NONchkvardims:
1536            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1537            varsadd.append(nonvd)
1538        ncvar.insert_variable(ncobj, 'clivi', diagout, diagoutd, diagoutvd, newnc)
1539
1540# WRFclwvi WRF water cloud-condensed path WRFdens, QCLOUD, QICE, QGRAUPEL, QHAIL
1541    elif diagn == 'WRFclwvi':
1542           
1543        var0 = WRFdens
1544        qtot = ncobj.variables[depvars[1]]
1545        qtotv = ncobj.variables[depvars[1]]
1546        Nspecies = len(depvars) - 2
1547        for iv in range(Nspecies):
1548            if ncobj.variables.has_key(depvars[iv+2]):
1549                var1 = ncobj.variables[depvars[iv+2]]
1550                qtotv = qtotv + var1[:]
1551
1552        dnamesvar = list(qtot.dimensions)
1553        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1554
1555        diagout, diagoutd, diagoutvd = diag.compute_clwvl(var0, qtotv, dnamesvar,dvnamesvar)
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        ncvar.insert_variable(ncobj, 'clwvi', diagout, diagoutd, diagoutvd, newnc)
1564
1565# WRF_denszint WRF vertical integration as WRFdens, sum(Q[water species1], ..., Q[water speciesN]), varn=[varN]
1566    elif diagn == 'WRF_denszint':
1567           
1568        var0 = WRFdens
1569        varn = depvars[1].split('=')[1]
1570        qtot = ncobj.variables[depvars[2]]
1571        qtotv = ncobj.variables[depvars[2]]
1572        Nspecies = len(depvars) - 2
1573        for iv in range(Nspecies):
1574            if ncobj.variables.has_key(depvars[iv+2]):
1575                var1 = ncobj.variables[depvars[iv+2]]
1576                qtotv = qtotv + var1[:]
1577
1578        dnamesvar = list(qtot.dimensions)
1579        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1580
1581        diagout, diagoutd, diagoutvd = diag.compute_clwvl(var0, qtotv, dnamesvar,dvnamesvar)
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        ncvar.insert_variable(ncobj, varn, diagout, diagoutd, diagoutvd, newnc)
1590
1591# WRFgeop geopotential from WRF as PH + PHB
1592    elif diagn == 'WRFgeop':
1593        var0 = ncobj.variables[depvars[0]][:]
1594        var1 = ncobj.variables[depvars[1]][:]
1595
1596        # de-staggering geopotential
1597        diagout0 = var0 + var1
1598        dt = diagout0.shape[0]
1599        dz = diagout0.shape[1]
1600        dy = diagout0.shape[2]
1601        dx = diagout0.shape[3]
1602
1603        diagout = np.zeros((dt,dz-1,dy,dx), dtype=np.float)
1604        diagout = 0.5*(diagout0[:,1:dz,:,:]+diagout0[:,0:dz-1,:,:])
1605
1606        # Removing the nonChecking variable-dimensions from the initial list
1607        varsadd = []
1608        diagoutvd = list(dvnames)
1609        for nonvd in NONchkvardims:
1610            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1611            varsadd.append(nonvd)
1612
1613        ncvar.insert_variable(ncobj, 'zg', diagout, dnames, diagoutvd, newnc)
1614
1615# WRFpotevap_orPM potential evapotranspiration following Penman-Monteith formulation
1616#   implemented in ORCHIDEE (in src_sechiba/enerbil.f90) as: WRFdens, UST, U10, V10, T2, PSFC, QVAPOR
1617    elif diagn == 'WRFpotevap_orPM':
1618        var0 = WRFdens[:,0,:,:]
1619        var1 = ncobj.variables[depvars[1]][:]
1620        var2 = ncobj.variables[depvars[2]][:]
1621        var3 = ncobj.variables[depvars[3]][:]
1622        var4 = ncobj.variables[depvars[4]][:]
1623        var5 = ncobj.variables[depvars[5]][:]
1624        var6 = ncobj.variables[depvars[6]][:,0,:,:]
1625
1626        dnamesvar = list(ncobj.variables[depvars[1]].dimensions)
1627        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1628
1629        diagout = np.zeros(var1.shape, dtype=np.float)
1630        diagout, diagoutd, diagoutvd = diag.Forcompute_potevap_orPM(var0, var1, var2,\
1631          var3, var4, var5, var6, dnamesvar, dvnamesvar)
1632
1633        # Removing the nonChecking variable-dimensions from the initial list
1634        varsadd = []
1635        for nonvd in NONchkvardims:
1636            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1637            varsadd.append(nonvd)
1638
1639        ncvar.insert_variable(ncobj, 'evspsblpot', diagout, diagoutd, diagoutvd, newnc)
1640
1641# WRFmslp_ptarget sea-level pressure following ECMWF method as PSFC, HGT, WRFt, WRFp, ZNU, ZNW
1642    elif diagn == 'WRFpsl_ecmwf':
1643        var0 = ncobj.variables[depvars[0]][:]
1644        var1 = ncobj.variables[depvars[1]][0,:,:]
1645        var2 = WRFt[:,0,:,:]
1646        var4 = WRFp[:,0,:,:]
1647        var5 = ncobj.variables[depvars[4]][0,:]
1648        var6 = ncobj.variables[depvars[5]][0,:]
1649
1650        # This is quite too appriximate!! passing pressure at half-levels to 2nd full
1651        #   level, using eta values at full (ZNW) and half (ZNU) mass levels
1652        var3 = WRFp[:,0,:,:] + (var6[1] - var5[0])*(WRFp[:,1,:,:] - WRFp[:,0,:,:])/  \
1653          (var5[1]-var5[0])
1654
1655        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1656        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1657
1658        diagout = np.zeros(var0.shape, dtype=np.float)
1659        diagout, diagoutd, diagoutvd = diag.Forcompute_psl_ecmwf(var0, var1, var2,   \
1660          var3, var4, dnamesvar, dvnamesvar)
1661
1662        # Removing the nonChecking variable-dimensions from the initial list
1663        varsadd = []
1664        for nonvd in NONchkvardims:
1665            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1666            varsadd.append(nonvd)
1667
1668        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
1669
1670# WRFmslp_ptarget sea-level pressure following ptarget method as WRFp, PSFC, WRFt, HGT, QVAPOR
1671    elif diagn == 'WRFpsl_ptarget':
1672        var0 = WRFp
1673        var1 = ncobj.variables[depvars[1]][:]
1674        var2 = WRFt
1675        var3 = ncobj.variables[depvars[3]][0,:,:]
1676        var4 = ncobj.variables[depvars[4]][:]
1677
1678        dnamesvar = list(ncobj.variables[depvars[4]].dimensions)
1679        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1680
1681        diagout = np.zeros(var0.shape, dtype=np.float)
1682        diagout, diagoutd, diagoutvd = diag.Forcompute_psl_ptarget(var0, var1, var2, \
1683          var3, var4, 700000., dnamesvar, dvnamesvar)
1684
1685        # Removing the nonChecking variable-dimensions from the initial list
1686        varsadd = []
1687        for nonvd in NONchkvardims:
1688            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1689            varsadd.append(nonvd)
1690
1691        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
1692
1693# WRFp pressure from WRF as P + PB
1694    elif diagn == 'WRFp':
1695        var0 = ncobj.variables[depvars[0]][:]
1696        var1 = ncobj.variables[depvars[1]][:]
1697           
1698        diagout = var0 + var1
1699        diagoutd = list(ncobj.variables[depvars[0]].dimensions)
1700        diagoutvd = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1701
1702        # Removing the nonChecking variable-dimensions from the initial list
1703        varsadd = []
1704        for nonvd in NONchkvardims:
1705            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1706            varsadd.append(nonvd)
1707
1708        ncvar.insert_variable(ncobj, 'pres', diagout, diagoutd, diagoutvd, newnc)
1709
1710# WRFpos
1711    elif diagn == 'WRFpos':
1712           
1713        dnamesvar = ncobj.variables['MAPFAC_M'].dimensions
1714        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1715
1716        ncvar.insert_variable(ncobj, 'WRFpos', WRFpos, dnamesvar, dvnamesvar, newnc)
1717
1718# WRFprw WRF water vapour path WRFdens, QVAPOR
1719    elif diagn == 'WRFprw':
1720           
1721        var0 = WRFdens
1722        var1 = ncobj.variables[depvars[1]]
1723
1724        dnamesvar = list(var1.dimensions)
1725        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1726
1727        diagout, diagoutd, diagoutvd = diag.compute_prw(var0, var1, dnamesvar,dvnamesvar)
1728
1729        # Removing the nonChecking variable-dimensions from the initial list
1730        varsadd = []
1731        diagoutvd = list(dvnames)
1732        for nonvd in NONchkvardims:
1733            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1734            varsadd.append(nonvd)
1735        ncvar.insert_variable(ncobj, 'prw', diagout, diagoutd, diagoutvd, newnc)
1736
1737# WRFrh (P, T, QVAPOR)
1738    elif diagn == 'WRFrh':
1739           
1740        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
1741        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1742
1743        ncvar.insert_variable(ncobj, 'hur', WRFrh, dnames, dvnames, newnc)
1744
1745# WRFrhs (PSFC, T2, Q2)
1746    elif diagn == 'WRFrhs':
1747           
1748        var0 = ncobj.variables[depvars[0]][:]
1749        var1 = ncobj.variables[depvars[1]][:]
1750        var2 = ncobj.variables[depvars[2]][:]
1751
1752        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
1753        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1754
1755        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
1756        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
1757
1758# rvors (u10, v10, WRFpos)
1759    elif diagn == 'WRFrvors':
1760           
1761        var0 = ncobj.variables[depvars[0]]
1762        var1 = ncobj.variables[depvars[1]]
1763
1764        diagout = rotational_z(var0, var1, distx)
1765
1766        dnamesvar = ncobj.variables[depvars[0]].dimensions
1767        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1768
1769        ncvar.insert_variable(ncobj, 'rvors', diagout, dnamesvar, dvnamesvar, newnc)
1770
1771# WRFt (T, P, PB)
1772    elif diagn == 'WRFt':
1773        var0 = ncobj.variables[depvars[0]][:]
1774        var1 = ncobj.variables[depvars[1]][:]
1775        var2 = ncobj.variables[depvars[2]][:]
1776
1777        p0=100000.
1778        p=var1 + var2
1779
1780        WRFt = (var0 + 300.)*(p/p0)**(2./7.)
1781
1782        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1783        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1784
1785        # Removing the nonChecking variable-dimensions from the initial list
1786        varsadd = []
1787        diagoutvd = list(dvnames)
1788        for nonvd in NONchkvardims:
1789            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1790            varsadd.append(nonvd)
1791
1792        ncvar.insert_variable(ncobj, 'ta', WRFt, dnames, diagoutvd, newnc)
1793
1794# WRFtda (WRFrh, WRFt)
1795    elif diagn == 'WRFtda':
1796        ARM2 = fdef.module_definitions.arm2
1797        ARM3 = fdef.module_definitions.arm3
1798
1799        gammatarh = np.log(WRFrh) + ARM2*(WRFt-273.15)/((WRFt-273.15)+ARM3)
1800        td = ARM3*gammatarh/(ARM2-gammatarh)
1801
1802        dnamesvar = list(ncobj.variables['T'].dimensions)
1803        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1804
1805        # Removing the nonChecking variable-dimensions from the initial list
1806        varsadd = []
1807        diagoutvd = list(dvnames)
1808        for nonvd in NONchkvardims:
1809            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1810            varsadd.append(nonvd)
1811
1812        ncvar.insert_variable(ncobj, 'tda', td, dnames, diagoutvd, newnc)
1813
1814# WRFtdas (PSFC, T2, Q2)
1815    elif diagn == 'WRFtdas':
1816        ARM2 = fdef.module_definitions.arm2
1817        ARM3 = fdef.module_definitions.arm3
1818
1819        var0 = ncobj.variables[depvars[0]][:]
1820        var1 = ncobj.variables[depvars[1]][:]
1821        var2 = ncobj.variables[depvars[2]][:]
1822
1823        dnamesvar = list(ncobj.variables[depvars[1]].dimensions)
1824        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1825
1826        rhs, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
1827
1828        gammatarhs = np.log(rhs) + ARM2*(var1-273.15)/((var1-273.15)+ARM3)
1829        tdas = ARM3*gammatarhs/(ARM2-gammatarhs) + 273.15
1830
1831        # Removing the nonChecking variable-dimensions from the initial list
1832        varsadd = []
1833        diagoutvd = list(dvnames)
1834        for nonvd in NONchkvardims:
1835            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1836            varsadd.append(nonvd)
1837
1838        ncvar.insert_variable(ncobj, 'tdas', tdas, dnames, diagoutvd, newnc)
1839
1840# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
1841    elif diagn == 'WRFua':
1842        var0 = ncobj.variables[depvars[0]][:]
1843        var1 = ncobj.variables[depvars[1]][:]
1844        var2 = ncobj.variables[depvars[2]][:]
1845        var3 = ncobj.variables[depvars[3]][:]
1846
1847        # un-staggering variables
1848        if len(var0.shape) == 4:
1849            unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1850        elif len(var0.shape) == 3:
1851            # Asuming sunding point (dimt, dimz, dimstgx)
1852            unstgdims = [var0.shape[0], var0.shape[1]]
1853
1854        ua = np.zeros(tuple(unstgdims), dtype=np.float)
1855        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1856        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1857
1858        if len(var0.shape) == 4:
1859            unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1860            unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1861
1862            for iz in range(var0.shape[1]):
1863                ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
1864
1865            dnamesvar = ['Time','bottom_top','south_north','west_east']
1866
1867        elif len(var0.shape) == 3:
1868            unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1])
1869            unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1])
1870            for iz in range(var0.shape[1]):
1871                ua[:,iz] = unstgvar0[:,iz]*var3 - unstgvar1[:,iz]*var2
1872
1873            dnamesvar = ['Time','bottom_top']
1874
1875        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1876
1877        # Removing the nonChecking variable-dimensions from the initial list
1878        varsadd = []
1879        diagoutvd = list(dvnames)
1880        for nonvd in NONchkvardims:
1881            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1882            varsadd.append(nonvd)
1883
1884        ncvar.insert_variable(ncobj, 'ua', ua, dnames, diagoutvd, newnc)
1885
1886# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
1887    elif diagn == 'WRFva':
1888        var0 = ncobj.variables[depvars[0]][:]
1889        var1 = ncobj.variables[depvars[1]][:]
1890        var2 = ncobj.variables[depvars[2]][:]
1891        var3 = ncobj.variables[depvars[3]][:]
1892
1893        # un-staggering variables
1894        if len(var0.shape) == 4:
1895            unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1896        elif len(var0.shape) == 3:
1897            # Asuming sunding point (dimt, dimz, dimstgx)
1898            unstgdims = [var0.shape[0], var0.shape[1]]
1899
1900        va = np.zeros(tuple(unstgdims), dtype=np.float)
1901        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1902        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1903
1904        if len(var0.shape) == 4:
1905            unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1906            unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1907
1908            for iz in range(var0.shape[1]):
1909                va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
1910
1911            dnamesvar = ['Time','bottom_top','south_north','west_east']
1912
1913        elif len(var0.shape) == 3:
1914            unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1])
1915            unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1])
1916            for iz in range(var0.shape[1]):
1917                va[:,iz] = unstgvar0[:,iz]*var2 + unstgvar1[:,iz]*var3
1918
1919            dnamesvar = ['Time','bottom_top']
1920
1921        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1922
1923        # Removing the nonChecking variable-dimensions from the initial list
1924        varsadd = []
1925        diagoutvd = list(dvnames)
1926        for nonvd in NONchkvardims:
1927            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1928            varsadd.append(nonvd)
1929        ncvar.insert_variable(ncobj, 'va', va, dnames, diagoutvd, newnc)
1930
1931
1932# WRFwd (U, V, SINALPHA, COSALPHA) to be rotated !!
1933    elif diagn == 'WRFwd':
1934        var0 = ncobj.variables[depvars[0]][:]
1935        var1 = ncobj.variables[depvars[1]][:]
1936        var2 = ncobj.variables[depvars[2]][:]
1937        var3 = ncobj.variables[depvars[3]][:]
1938
1939        # un-staggering variables
1940        if len(var0.shape) == 4:
1941            unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1942        elif len(var0.shape) == 3:
1943            # Asuming sunding point (dimt, dimz, dimstgx)
1944            unstgdims = [var0.shape[0], var0.shape[1]]
1945
1946        ua = np.zeros(tuple(unstgdims), dtype=np.float)
1947        va = np.zeros(tuple(unstgdims), dtype=np.float)
1948        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1949        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1950
1951        if len(var0.shape) == 4:
1952            unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1953            unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1954
1955            for iz in range(var0.shape[1]):
1956                ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
1957                va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
1958
1959            dnamesvar = ['Time','bottom_top','south_north','west_east']
1960
1961        elif len(var0.shape) == 3:
1962            unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1])
1963            unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1])
1964            for iz in range(var0.shape[1]):
1965                ua[:,iz] = unstgvar0[:,iz]*var3 - unstgvar1[:,iz]*var2
1966                va[:,iz] = unstgvar0[:,iz]*var2 + unstgvar1[:,iz]*var3
1967
1968            dnamesvar = ['Time','bottom_top']
1969
1970        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1971
1972        wd, dnames, dvnames = diag.compute_wd(ua, va, dnamesvar, dvnamesvar)
1973
1974        # Removing the nonChecking variable-dimensions from the initial list
1975        varsadd = []
1976        diagoutvd = list(dvnames)
1977        for nonvd in NONchkvardims:
1978            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1979            varsadd.append(nonvd)
1980
1981        ncvar.insert_variable(ncobj, 'wd', wd, dnames, diagoutvd, newnc)
1982
1983# WRFtime
1984    elif diagn == 'WRFtime':
1985           
1986        diagout = WRFtime
1987
1988        dnamesvar = ['Time']
1989        dvnamesvar = ['Times']
1990
1991        ncvar.insert_variable(ncobj, 'time', diagout, dnamesvar, dvnamesvar, newnc)
1992
1993# ws (U, V)
1994    elif diagn == 'ws':
1995           
1996        # un-staggering variables
1997        if len(var0.shape) == 4:
1998            unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1999        elif len(var0.shape) == 3:
2000            # Asuming sunding point (dimt, dimz, dimstgx)
2001            unstgdims = [var0.shape[0], var0.shape[1]]
2002
2003        ua = np.zeros(tuple(unstgdims), dtype=np.float)
2004        va = np.zeros(tuple(unstgdims), dtype=np.float)
2005        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
2006        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
2007
2008        if len(var0.shape) == 4:
2009            unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
2010            unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
2011
2012            for iz in range(var0.shape[1]):
2013                ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
2014                va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
2015
2016            dnamesvar = ['Time','bottom_top','south_north','west_east']
2017
2018        elif len(var0.shape) == 3:
2019            unstgvar0 = 0.5*(var0[:,:,0] + var0[:,:,1])
2020            unstgvar1 = 0.5*(var1[:,:,0] + var1[:,:,1])
2021            for iz in range(var0.shape[1]):
2022                ua[:,iz] = unstgvar0[:,iz]*var3 - unstgvar1[:,iz]*var2
2023                va[:,iz] = unstgvar0[:,iz]*var2 + unstgvar1[:,iz]*var3
2024
2025            dnamesvar = ['Time','bottom_top']
2026
2027        diagout = np.sqrt(unstgvar0*unstgvar0 + unstgvar1*unstgvar1)
2028
2029#        dnamesvar = ncobj.variables[depvars[0]].dimensions
2030        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2031
2032        # Removing the nonChecking variable-dimensions from the initial list
2033        varsadd = []
2034        diagoutvd = list(dvnamesvar)
2035        for nonvd in NONchkvardims:
2036            if gen.searchInlist(dvnamesvar,nonvd): diagoutvd.remove(nonvd)
2037            varsadd.append(nonvd)
2038        ncvar.insert_variable(ncobj, 'ws', diagout, dnamesvar, diagoutvd, newnc)
2039
2040# wss (u10, v10)
2041    elif diagn == 'wss':
2042           
2043        var0 = ncobj.variables[depvars[0]][:]
2044        var1 = ncobj.variables[depvars[1]][:]
2045
2046        diagout = np.sqrt(var0*var0 + var1*var1)
2047
2048        dnamesvar = ncobj.variables[depvars[0]].dimensions
2049        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
2050
2051        ncvar.insert_variable(ncobj, 'wss', diagout, dnamesvar, dvnamesvar, newnc)
2052
2053# WRFheight height from WRF geopotential as WRFGeop/g
2054    elif diagn == 'WRFheight':
2055           
2056        diagout = WRFgeop/grav
2057
2058        # Removing the nonChecking variable-dimensions from the initial list
2059        varsadd = []
2060        diagoutvd = list(dvnames)
2061        for nonvd in NONchkvardims:
2062            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2063            varsadd.append(nonvd)
2064
2065        ncvar.insert_variable(ncobj, 'zhgt', diagout, dnames, diagoutvd, newnc)
2066
2067# WRFheightrel relative-height from WRF geopotential as WRFgeop(PH + PHB)/g-HGT 'WRFheightrel|PH@PHB@HGT
2068    elif diagn == 'WRFheightrel':
2069        var0 = ncobj.variables[depvars[0]][:]
2070        var1 = ncobj.variables[depvars[1]][:]
2071        var2 = ncobj.variables[depvars[2]][:]
2072
2073        dimz = var0.shape[1]
2074        diagout = np.zeros(tuple(var0.shape), dtype=np.float)
2075        for iz in range(dimz):
2076            diagout[:,iz,:,:] = (var0[:,iz,:,:]+ var1[:,iz,:,:])/grav - var2
2077
2078        # Removing the nonChecking variable-dimensions from the initial list
2079        varsadd = []
2080        diagoutvd = list(dvnames)
2081        for nonvd in NONchkvardims:
2082            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2083            varsadd.append(nonvd)
2084
2085        ncvar.insert_variable(ncobj, 'zhgtrel', diagout, dnames, diagoutvd, newnc)
2086
2087# WRFzmla_gen generic boundary layer hieght computation from WRF theta, QVAPOR, WRFgeop, HGT,
2088    elif diagn == 'WRFzmlagen':
2089        var0 = ncobj.variables[depvars[0]][:]+300.
2090        var1 = ncobj.variables[depvars[1]][:]
2091        dimz = var0.shape[1]
2092        var2 = WRFgeop[:,1:dimz+1,:,:]/9.8
2093        var3 = ncobj.variables[depvars[3]][0,:,:]
2094
2095        diagout, diagoutd, diagoutvd = diag.Forcompute_zmla_gen(var0,var1,var2,var3, \
2096          dnames,dvnames)
2097
2098        # Removing the nonChecking variable-dimensions from the initial list
2099        varsadd = []
2100        for nonvd in NONchkvardims:
2101            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2102            varsadd.append(nonvd)
2103
2104        ncvar.insert_variable(ncobj, 'zmla', diagout, diagoutd, diagoutvd, newnc)
2105
2106# WRFzwind wind extrapolation at a given height using power law computation from WRF
2107#   U, V, WRFz, U10, V10, SINALPHA, COSALPHA, z=[zval]
2108    elif diagn == 'WRFzwind':
2109        var0 = ncobj.variables[depvars[0]][:]
2110        var1 = ncobj.variables[depvars[1]][:]
2111        var2 = WRFz
2112        var3 = ncobj.variables[depvars[3]][:]
2113        var4 = ncobj.variables[depvars[4]][:]
2114        var5 = ncobj.variables[depvars[5]][0,:,:]
2115        var6 = ncobj.variables[depvars[6]][0,:,:]
2116        var7 = np.float(depvars[7].split('=')[1])
2117
2118        # un-staggering 3D winds
2119        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
2120        va = np.zeros(tuple(unstgdims), dtype=np.float)
2121        unvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
2122        unvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
2123        unvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
2124        unvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
2125
2126        diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwind(unvar0,      \
2127          unvar1, var2, var3, var4, var5, var6, var7, dnames, dvnames)
2128
2129        # Removing the nonChecking variable-dimensions from the initial list
2130        varsadd = []
2131        for nonvd in NONchkvardims:
2132            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2133            varsadd.append(nonvd)
2134
2135        ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
2136        ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)
2137
2138# WRFzwind wind extrapolation at a given hieght using logarithmic law computation
2139#   from WRF U, V, WRFz, U10, V10, SINALPHA, COSALPHA, z=[zval]
2140    elif diagn == 'WRFzwind_log':
2141        var0 = ncobj.variables[depvars[0]][:]
2142        var1 = ncobj.variables[depvars[1]][:]
2143        var2 = WRFz
2144        var3 = ncobj.variables[depvars[3]][:]
2145        var4 = ncobj.variables[depvars[4]][:]
2146        var5 = ncobj.variables[depvars[5]][0,:,:]
2147        var6 = ncobj.variables[depvars[6]][0,:,:]
2148        var7 = np.float(depvars[7].split('=')[1])
2149
2150        # un-staggering 3D winds
2151        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
2152        va = np.zeros(tuple(unstgdims), dtype=np.float)
2153        unvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
2154        unvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
2155        unvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
2156        unvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
2157
2158        diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwind_log(unvar0,  \
2159          unvar1, var2, var3, var4, var5, var6, var7, dnames, dvnames)
2160
2161        # Removing the nonChecking variable-dimensions from the initial list
2162        varsadd = []
2163        for nonvd in NONchkvardims:
2164            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2165            varsadd.append(nonvd)
2166
2167        ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
2168        ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)
2169
2170# WRFzwindMO wind extrapolation at a given height computation using Monin-Obukhow
2171#   theory from WRF UST, ZNT, RMOL, U10, V10, SINALPHA, COSALPHA, z=[zval]
2172#   NOTE: only useful for [zval] < 80. m
2173    elif diagn == 'WRFzwindMO':
2174        var0 = ncobj.variables[depvars[0]][:]
2175        var1 = ncobj.variables[depvars[1]][:]
2176        var2 = ncobj.variables[depvars[2]][:]
2177        var3 = ncobj.variables[depvars[3]][:]
2178        var4 = ncobj.variables[depvars[4]][:]
2179        var5 = ncobj.variables[depvars[5]][0,:,:]
2180        var6 = ncobj.variables[depvars[6]][0,:,:]
2181        var7 = np.float(depvars[7].split('=')[1])
2182
2183        diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwindMO(var0, var1,\
2184          var2, var3, var4, var5, var6, var7, dnames, dvnames)
2185
2186        # Removing the nonChecking variable-dimensions from the initial list
2187        varsadd = []
2188        for nonvd in NONchkvardims:
2189            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2190            varsadd.append(nonvd)
2191
2192        ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
2193        ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)
2194
2195# zmla_gen generic boundary layer hieght computation from generic 2D-file theta, qv, zg, orog,
2196    elif diagn == 'zmlagen':
2197        var0 = ncobj.variables[depvars[0]][:]
2198        var1 = ncobj.variables[depvars[1]][:]
2199        dimz = var0.shape[1]
2200        var2 = ncobj.variables[depvars[2]][:]
2201        var3 = ncobj.variables[depvars[3]][:]
2202
2203        diagout, diagoutd, diagoutvd = diag.Forcompute_zmla_gen(var0,var1,var2,var3, \
2204          dnames,dvnames)
2205
2206        # Removing the nonChecking variable-dimensions from the initial list
2207        varsadd = []
2208        for nonvd in NONchkvardims:
2209            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2210            varsadd.append(nonvd)
2211
2212        ncvar.insert_variable(ncobj, 'zmla', diagout, diagoutd, diagoutvd, newnc)
2213
2214# zmla_genUWsnd generic boundary layer hieght computation from UWyoming sounding file theta, qv, zg
2215    elif diagn == 'zmlagenUWsnd':
2216        var0 = ncobj.variables[depvars[0]][:]
2217        var1 = ncobj.variables[depvars[1]][:]
2218        dimz = var0.shape[1]
2219        var2 = ncobj.variables[depvars[2]][:]/9.8
2220        var3 = ncobj.getncattr('Station_elevation')
2221
2222        diagout, diagoutd, diagoutvd = diag.Forcompute_zmla_gen(var0,var1,var2,var3, \
2223          dnames,dvnames)
2224
2225        # No need to remove topography height
2226        diagout = diagout + var3
2227
2228        # Removing the nonChecking variable-dimensions from the initial list
2229        varsadd = []
2230        for nonvd in NONchkvardims:
2231            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
2232            varsadd.append(nonvd)
2233
2234        ncvar.insert_variable(ncobj, 'zmla', diagout, diagoutd, diagoutvd, newnc)
2235
2236    else:
2237        print errormsg
2238        print '  ' + main + ": diagnostic '" + diagn + "' not ready!!!"
2239        print '    available diagnostics: ', availdiags
2240        quit(-1)
2241
2242    newnc.sync()
2243    # Adding that additional variables required to compute some diagnostics which
2244    #   where not in the original file
2245    print '  adding additional variables...'
2246    for vadd in varsadd:
2247        if not gen.searchInlist(newnc.variables.keys(),vadd) and                     \
2248          dictcompvars.has_key(vadd):
2249            attrs = dictcompvars[vadd]
2250            vvn = attrs['name']
2251            if not gen.searchInlist(newnc.variables.keys(), vvn):
2252                iidvn = dvnames.index(vadd)
2253                dnn = dnames[iidvn]
2254                if vadd == 'WRFtime':
2255                    dvarvals = WRFtime[:]
2256                newvar = newnc.createVariable(vvn, 'f8', (dnn))
2257                newvar[:] = dvarvals
2258                for attn in attrs.keys():
2259                    if attn != 'name':
2260                        attv = attrs[attn]
2261                        ncvar.set_attribute(newvar, attn, attv)
2262
2263#   end of diagnostics
2264
2265# Global attributes
2266##
2267ncvar.add_global_PyNCplot(newnc, main, None, '2.0')
2268
2269gorigattrs = ncobj.ncattrs()
2270for attr in gorigattrs:
2271    attrv = ncobj.getncattr(attr)
2272    atvar = ncvar.set_attribute(newnc, attr, attrv)
2273
2274ncobj.close()
2275newnc.close()
2276
2277print '\n' + main + ': successfull writting of diagnostics file "' + ofile + '" !!!'
Note: See TracBrowser for help on using the repository browser.