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

Last change on this file since 2201 was 2141, checked in by lfita, 6 years ago

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