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

Last change on this file since 2680 was 2679, checked in by lfita, 6 years ago

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