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

Last change on this file since 2617 was 2390, checked in by lfita, 6 years ago

Adding:

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