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

Last change on this file since 1808 was 1806, checked in by lfita, 7 years ago

Adding `WRFpotevap_orPM' to the list of diagnostics
File size: 52.9 KB
Line 
1# Python script to comput diagnostics
2# L. Fita, LMD. CNR, UPMC-Jussieu, Paris, France
3# File diagnostics.inf provides the combination of variables to get the desired diagnostic
4#   To be used with module_ForDiagnostics.F90, module_ForDiagnosticsVars.F90, module_generic.F90
5#      foudre: f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90 >& run_f2py.log
6#      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
7
8## 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@XTIME' -f WRF_LMDZ/NPv31/wrfout_d01_1980-03-01_00:00:00
9## e.g. # diagnostics.py -f /home/lluis/PY/diagnostics.inf -d variable_combo -v WRFprc
10
11# Available general pupose diagnostics (model independent) providing (varv1, varv2, ..., dimns, dimvns)
12# compute_accum: Function to compute the accumulation of a variable
13# compute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following
14#   newmicro.F90 from LMDZ compute_clt(cldfra, pres, dimns, dimvns)
15# compute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ
16# compute_clivi: Function to compute cloud-ice water path (clivi)
17# compute_clwvl: Function to compute condensed water path (clwvl)
18# compute_deaccum: Function to compute the deaccumulation of a variable
19# compute_mslp: Function to compute mslp: mean sea level pressure following p_interp.F90 from WRF
20# compute_OMEGAw: Function to transform OMEGA [Pas-1] to velocities [ms-1]
21# compute_prw: Function to compute water vapour path (prw)
22# compute_rh: Function to compute relative humidity following 'Tetens' equation (T,P) ...'
23# compute_td: Function to compute the dew point temperature
24# compute_turbulence: Function to compute the rubulence term of the Taylor's decomposition ...'
25# compute_wds: Function to compute the wind direction
26# compute_wss: Function to compute the wind speed
27# compute_WRFuava: Function to compute geographical rotated WRF 3D winds
28# compute_WRFuasvas: Fucntion to compute geographical rotated WRF 2-meter winds
29# derivate_centered: Function to compute the centered derivate of a given field
30# def Forcompute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ via Fortran subroutine
31# Forcompute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ via a Fortran module
32# Forcompute_psl_ptarget: Function to compute the sea-level pressure following target_pressure value found in `p_interp.F'
33
34# Others just providing variable values
35# var_cllmh: Fcuntion to compute cllmh on a 1D column
36# var_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ using 1D vertical column values
37# var_mslp: Fcuntion to compute mean sea-level pressure
38# var_virtualTemp: This function returns virtual temperature in K,
39# var_WRFtime: Function to copmute CFtimes from WRFtime variable
40# rotational_z: z-component of the rotatinoal of horizontal vectorial field
41# turbulence_var: Function to compute the Taylor's decomposition turbulence term from a a given variable
42
43from optparse import OptionParser
44import numpy as np
45from netCDF4 import Dataset as NetCDFFile
46import os
47import re
48import nc_var_tools as ncvar
49import generic_tools as gen
50import datetime as dtime
51import module_ForDiag as fdin
52import diag_tools as diag
53
54main = 'diagnostics.py'
55errormsg = 'ERROR -- error -- ERROR -- error'
56warnmsg = 'WARNING -- warning -- WARNING -- warning'
57
58# Constants
59grav = 9.81
60
61
62####### ###### ##### #### ### ## #
63comboinf="\nIF -d 'variable_combo', provides information of the combination to obtain -v [varn] with the ASCII file with the combinations as -f [combofile]"
64
65parser = OptionParser()
66parser.add_option("-f", "--netCDF_file", dest="ncfile", help="file to use", metavar="FILE")
67parser.add_option("-d", "--dimensions", dest="dimns", 
68  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, 
69  metavar="LABELS")
70parser.add_option("-v", "--variables", dest="varns", 
71  help=" [varn1]|[var11]@[...[varN1]],[...,[varnM]|[var1M]@[...[varLM]]] ',' list of variables to compute [varnK] and its necessary ones [var1K]...[varPK]", metavar="VALUES")
72
73(opts, args) = parser.parse_args()
74
75#######    #######
76## MAIN
77    #######
78availdiags = ['ACRAINTOT', 'accum', 'clt', 'cllmh', 'deaccum', 'LMDZrh', 'mslp',     \
79  'OMEGAw', 'RAINTOT',                                                               \
80  'rvors', 'td', 'turbulence', 'WRFcape_afwa', 'WRFclivi', 'WRFclwvi', 'WRFgeop',    \
81  'WRFmrso', 'WRFpotevap_orPM', 'WRFp', 'WRFpsl_ecmwf',                              \
82  'WRFpsl_ptarget', 'WRFrvors', 'WRFslw', 'ws', 'wds', 'wss', 'WRFheight',           \
83  'WRFheightrel', 'WRFua', 'WRFva', 'WRFzwind', 'WRFzwind_log', 'WRFzwindMO']
84
85methods = ['accum', 'deaccum']
86
87# Variables not to check
88NONcheckingvars = ['cllmh', 'deaccum', 'TSrhs', 'TStd', 'TSwds', 'TSwss', 'WRFbils', \
89  'WRFclivi', 'WRFclwvi', 'WRFdens', 'WRFgeop',                                      \
90  'WRFp', 'WRFtd',                                                                   \
91  'WRFpos', 'WRFpotevap_orPM', 'WRFprc', 'WRFprls', 'WRFrh', 'LMDZrh', 'LMDZrhs',    \
92  'WRFrhs', 'WRFrvors',                                                              \
93  'WRFt', 'WRFtime', 'WRFua', 'WRFva', 'WRFwds', 'WRFwss', 'WRFheight', 'WRFz']
94
95NONchkvardims = ['WRFtime']
96
97ofile = 'diagnostics.nc'
98
99dimns = opts.dimns
100varns = opts.varns
101
102# Special method. knowing variable combination
103##
104if opts.dimns == 'variable_combo':
105    print warnmsg
106    print '  ' + main + ': knowing variable combination !!!'
107    combination = variable_combo(opts.varns,opts.ncfile)
108    print '     COMBO: ' + combination
109    quit(-1)
110
111if not os.path.isfile(opts.ncfile):
112    print errormsg
113    print '   ' + main + ": file '" + opts.ncfile + "' does not exist !!"
114    quit(-1)
115
116ncobj = NetCDFFile(opts.ncfile, 'r')
117
118# Looking for specific variables that might be use in more than one diagnostic
119WRFgeop_compute = False
120WRFp_compute = False
121WRFt_compute = False
122WRFrh_compute = False
123WRFght_compute = False
124WRFdens_compute = False
125WRFpos_compute = False
126WRFtime_compute = False
127WRFz_compute = False
128
129# File creation
130newnc = NetCDFFile(ofile,'w')
131
132# dimensions
133dimvalues = dimns.split(',')
134dnames = []
135dvnames = []
136
137for dimval in dimvalues:
138    dn = dimval.split('@')[0]
139    dnv = dimval.split('@')[1]
140    dnames.append(dn)
141    dvnames.append(dnv)
142    # Is there any dimension-variable which should be computed?
143    if dnv == 'WRFgeop':WRFgeop_compute = True
144    if dnv == 'WRFp': WRFp_compute = True
145    if dnv == 'WRFt': WRFt_compute = True
146    if dnv == 'WRFrh': WRFrh_compute = True
147    if dnv == 'WRFght': WRFght_compute = True
148    if dnv == 'WRFdens': WRFdens_compute = True
149    if dnv == 'WRFpos': WRFpos_compute = True
150    if dnv == 'WRFtime': WRFtime_compute = True
151    if dnv == 'WRFz':WRFz_compute = True
152
153# diagnostics to compute
154diags = varns.split(',')
155Ndiags = len(diags)
156
157for idiag in range(Ndiags):
158    if diags[idiag].split('|')[1].find('@') == -1:
159        depvars = diags[idiag].split('|')[1]
160        if depvars == 'WRFgeop':WRFgeop_compute = True
161        if depvars == 'WRFp': WRFp_compute = True
162        if depvars == 'WRFt': WRFt_compute = True
163        if depvars == 'WRFrh': WRFrh_compute = True
164        if depvars == 'WRFght': WRFght_compute = True
165        if depvars == 'WRFdens': WRFdens_compute = True
166        if depvars == 'WRFpos': WRFpos_compute = True
167        if depvars == 'WRFtime': WRFtime_compute = True
168        if depvars == 'WRFz': WRFz_compute = True
169    else:
170        depvars = diags[idiag].split('|')[1].split('@')
171        if gen.searchInlist(depvars, 'WRFgeop'): WRFgeop_compute = True
172        if gen.searchInlist(depvars, 'WRFp'): WRFp_compute = True
173        if gen.searchInlist(depvars, 'WRFt'): WRFt_compute = True
174        if gen.searchInlist(depvars, 'WRFrh'): WRFrh_compute = True
175        if gen.searchInlist(depvars, 'WRFght'): WRFght_compute = True
176        if gen.searchInlist(depvars, 'WRFdens'): WRFdens_compute = True
177        if gen.searchInlist(depvars, 'WRFpos'): WRFpos_compute = True
178        if gen.searchInlist(depvars, 'WRFtime'): WRFtime_compute = True
179        if gen.searchInlist(depvars, 'WRFz'): WRFz_compute = True
180
181# Dictionary with the new computed variables to be able to add them
182dictcompvars = {}
183if WRFgeop_compute:
184    print '  ' + main + ': Retrieving geopotential value from WRF as PH + PHB'
185    dimv = ncobj.variables['PH'].shape
186    WRFgeop = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
187
188    # Attributes of the variable
189    Vvals = gen.variables_values('WRFgeop')
190    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
191      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
192
193if WRFp_compute:
194    print '  ' + main + ': Retrieving pressure value from WRF as P + PB'
195    dimv = ncobj.variables['P'].shape
196    WRFp = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
197
198    # Attributes of the variable
199    Vvals = gen.variables_values('WRFp')
200    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
201      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
202
203if WRFght_compute:
204    print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...' 
205    WRFght = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
206
207    # Attributes of the variable
208    Vvals = gen.variables_values('WRFght')
209    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
210      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
211
212if WRFrh_compute:
213    print '    ' + main + ": computing relative humidity from WRF as 'Tetens'" +     \
214      ' equation (T,P) ...'
215    p0=100000.
216    p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
217    tk = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
218    qv = ncobj.variables['QVAPOR'][:]
219
220    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
221    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
222
223    WRFrh = qv/data2
224
225    # Attributes of the variable
226    Vvals = gen.variables_values('WRFrh')
227    dictcompvars['WRFrh'] = {'name': Vvals[0], 'standard_name': Vvals[1],            \
228      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
229
230if WRFt_compute:
231    print '    ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
232    p0=100000.
233    p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
234
235    WRFt = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
236
237    # Attributes of the variable
238    Vvals = gen.variables_values('WRFt')
239    dictcompvars['WRFt'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
240      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
241
242if WRFdens_compute:
243    print '    ' + main + ': computing air density from WRF as ((MU + MUB) * ' +     \
244      'DNW)/g ...'
245
246# Just we need in in absolute values: Size of the central grid cell
247##    dxval = ncobj.getncattr('DX')
248##    dyval = ncobj.getncattr('DY')
249##    mapfac = ncobj.variables['MAPFAC_M'][:]
250##    area = dxval*dyval*mapfac
251
252    mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
253    dnw = ncobj.variables['DNW'][:]
254
255    WRFdens = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]),        \
256      dtype=np.float)
257    levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)
258
259    for it in range(mu.shape[0]):
260        for iz in range(dnw.shape[1]):
261            levval.fill(np.abs(dnw[it,iz]))
262            WRFdens[it,iz,:,:] = levval
263            WRFdens[it,iz,:,:] = mu[it,:,:]*WRFdens[it,iz,:,:]/grav
264
265    # Attributes of the variable
266    Vvals = gen.variables_values('WRFdens')
267    dictcompvars['WRFdens'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
268      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
269
270if WRFpos_compute:
271# WRF positions from the lowest-leftest corner of the matrix
272    print '    ' + main + ': computing position from MAPFAC_M as sqrt(DY*j**2 + ' +  \
273      'DX*x**2)*MAPFAC_M ...'
274
275    mapfac = ncobj.variables['MAPFAC_M'][:]
276
277    distx = np.float(ncobj.getncattr('DX'))
278    disty = np.float(ncobj.getncattr('DY'))
279
280    print 'distx:',distx,'disty:',disty
281
282    dx = mapfac.shape[2]
283    dy = mapfac.shape[1]
284    dt = mapfac.shape[0]
285
286    WRFpos = np.zeros((dt, dy, dx), dtype=np.float)
287
288    for i in range(1,dx):
289        WRFpos[0,0,i] = distx*i/mapfac[0,0,i]
290    for j in range(1,dy):
291        i=0
292        WRFpos[0,j,i] = WRFpos[0,j-1,i] + disty/mapfac[0,j,i]
293        for i in range(1,dx):
294#            WRFpos[0,j,i] = np.sqrt((disty*j)**2. + (distx*i)**2.)/mapfac[0,j,i]
295#            WRFpos[0,j,i] = np.sqrt((disty*j)**2. + (distx*i)**2.)
296             WRFpos[0,j,i] = WRFpos[0,j,i-1] + distx/mapfac[0,j,i]
297
298    for it in range(1,dt):
299        WRFpos[it,:,:] = WRFpos[0,:,:]
300
301if WRFtime_compute:
302    print '    ' + main + ': computing time from WRF as CFtime(Times) ...'
303
304    refdate='19491201000000'
305    tunitsval='minutes'
306
307    timeobj = ncobj.variables['Times']
308    timewrfv = timeobj[:]
309
310    yrref=refdate[0:4]
311    monref=refdate[4:6]
312    dayref=refdate[6:8]
313    horref=refdate[8:10]
314    minref=refdate[10:12]
315    secref=refdate[12:14]
316
317    refdateS = yrref + '-' + monref + '-' + dayref + ' ' + horref + ':' + minref +   \
318      ':' + secref
319
320    dt = timeobj.shape[0]
321    WRFtime = np.zeros((dt), dtype=np.float)
322
323    for it in range(dt):
324        wrfdates = gen.datetimeStr_conversion(timewrfv[it,:],'WRFdatetime', 'matYmdHMS')
325        WRFtime[it] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
326
327    tunits = tunitsval + ' since ' + refdateS
328
329    # Attributes of the variable
330    dictcompvars['WRFtime'] = {'name': 'time', 'standard_name': 'time',              \
331      'long_name': 'time', 'units': tunits, 'calendar': 'gregorian'}
332
333if WRFz_compute:
334    print '  ' + main + ': Retrieving z: height above surface value from WRF as ' +  \
335      'unstagger(PH + PHB)/9.8-hgt'
336    dimv = ncobj.variables['PH'].shape
337    WRFzg = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/9.8
338
339    unzgd = (dimv[0], dimv[1]-1, dimv[2], dimv[3])
340    unzg = np.zeros(unzgd, dtype=np.float)
341    unzg = 0.5*(WRFzg[:,0:dimv[1]-1,:,:] + WRFzg[:,1:dimv[1],:,:])
342
343    WRFz = np.zeros(unzgd, dtype=np.float)
344    for iz in range(dimv[1]-1):
345        WRFz[:,iz,:,:] = unzg[:,iz,:,:] - ncobj.variables['HGT'][:]
346
347    # Attributes of the variable
348    Vvals = gen.variables_values('WRFz')
349    dictcompvars['WRFz'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
350      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
351
352### ## #
353# Going for the diagnostics
354### ## #
355print '  ' + main + ' ...'
356varsadd = []
357
358for idiag in range(Ndiags):
359    print '    diagnostic:',diags[idiag]
360    diagn = diags[idiag].split('|')[0]
361    depvars = diags[idiag].split('|')[1].split('@')
362    if diags[idiag].split('|')[1].find('@') != -1:
363        depvars = diags[idiag].split('|')[1].split('@')
364        if depvars[0] == 'deaccum': diagn='deaccum'
365        if depvars[0] == 'accum': diagn='accum'
366        for depv in depvars:
367            if not ncobj.variables.has_key(depv) and not                             \
368              gen.searchInlist(NONcheckingvars, depv) and                            \
369              not gen.searchInlist(methods, depv) and not depvars[0] == 'deaccum'    \
370              and not depvars[0] == 'accum' and not depv[0:2] == 'z=':
371                print errormsg
372                print '  ' + main + ": file '" + opts.ncfile +                       \
373                  "' does not have variable '" + depv + "' !!"
374                quit(-1)
375    else:
376        depvars = diags[idiag].split('|')[1]
377        if not ncobj.variables.has_key(depvars) and not                              \
378          gen.searchInlist(NONcheckingvars, depvars) and                             \
379          not gen.searchInlist(methods, depvars):
380            print errormsg
381            print '  ' + main + ": file '" + opts.ncfile +                           \
382              "' does not have variable '" + depvars + "' !!"
383            quit(-1)
384
385    print "\n    Computing '" + diagn + "' from: ", depvars, '...'
386
387# acraintot: accumulated total precipitation from WRF RAINC, RAINNC
388    if diagn == 'ACRAINTOT':
389           
390        var0 = ncobj.variables[depvars[0]]
391        var1 = ncobj.variables[depvars[1]]
392        diagout = var0[:] + var1[:]
393
394        dnamesvar = var0.dimensions
395        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
396
397        # Removing the nonChecking variable-dimensions from the initial list
398        varsadd = []
399        for nonvd in NONchkvardims:
400            if gen.searchInlist(dvnamesvar,nonvd): dvnamesvar.remove(nonvd)
401            varsadd.append(nonvd)
402
403        ncvar.insert_variable(ncobj, 'pracc', diagout, dnamesvar, dvnamesvar, newnc)
404
405# accum: acumulation of any variable as (Variable, time [as [tunits]
406#   from/since ....], newvarname)
407    elif diagn == 'accum':
408
409        var0 = ncobj.variables[depvars[0]]
410        var1 = ncobj.variables[depvars[1]]
411
412        dnamesvar = var0.dimensions
413        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
414
415        diagout, diagoutd, diagoutvd = diag.compute_accum(var0,dnamesvar,dvnamesvar)
416
417        CFvarn = ncvar.variables_values(depvars[0])[0]
418
419# Removing the flux
420        if depvars[1] == 'XTIME':
421            dtimeunits = var1.getncattr('description')
422            tunits = dtimeunits.split(' ')[0]
423        else:
424            dtimeunits = var1.getncattr('units')
425            tunits = dtimeunits.split(' ')[0]
426
427        dtime = (var1[1] - var1[0])*timeunits_seconds(tunits)
428
429        ncvar.insert_variable(ncobj, CFvarn + 'acc', diagout*dtime, diagoutd, diagoutvd, newnc)
430
431# cllmh with cldfra, pres
432    elif diagn == 'cllmh':
433           
434        var0 = ncobj.variables[depvars[0]]
435        if depvars[1] == 'WRFp':
436            var1 = WRFp
437        else:
438            var01 = ncobj.variables[depvars[1]]
439            if len(size(var1.shape)) < len(size(var0.shape)):
440                var1 = np.brodcast_arrays(var01,var0)[0]
441            else:
442                var1 = var01
443
444        diagout, diagoutd, diagoutvd = diag.Forcompute_cllmh(var0,var1,dnames,dvnames)
445
446        # Removing the nonChecking variable-dimensions from the initial list
447        varsadd = []
448        for nonvd in NONchkvardims:
449            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
450            varsadd.append(nonvd)
451
452        ncvar.insert_variable(ncobj, 'cll', diagout[0,:], diagoutd, diagoutvd, newnc)
453        ncvar.insert_variable(ncobj, 'clm', diagout[1,:], diagoutd, diagoutvd, newnc)
454        ncvar.insert_variable(ncobj, 'clh', diagout[2,:], diagoutd, diagoutvd, newnc)
455
456# clt with cldfra
457    elif diagn == 'clt':
458           
459        var0 = ncobj.variables[depvars]
460        diagout, diagoutd, diagoutvd = diag.Forcompute_clt(var0,dnames,dvnames)
461
462        # Removing the nonChecking variable-dimensions from the initial list
463        varsadd = []
464        for nonvd in NONchkvardims:
465            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
466            varsadd.append(nonvd)
467           
468        ncvar.insert_variable(ncobj, 'clt', diagout, diagoutd, diagoutvd, newnc)
469
470# deaccum: deacumulation of any variable as (Variable, time [as [tunits]
471#   from/since ....], newvarname)
472    elif diagn == 'deaccum':
473
474        var0 = ncobj.variables[depvars[1]]
475        var1 = ncobj.variables[depvars[2]]
476
477        dnamesvar = var0.dimensions
478        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
479
480        diagout, diagoutd, diagoutvd = diag.compute_deaccum(var0,dnamesvar,dvnamesvar)
481
482# Transforming to a flux
483        if depvars[2] == 'XTIME':
484            dtimeunits = var1.getncattr('description')
485            tunits = dtimeunits.split(' ')[0]
486        else:
487            dtimeunits = var1.getncattr('units')
488            tunits = dtimeunits.split(' ')[0]
489
490        dtime = (var1[1] - var1[0])*timeunits_seconds(tunits)
491        ncvar.insert_variable(ncobj, depvars[3], diagout/dtime, diagoutd, diagoutvd, newnc)
492
493# LMDZrh (pres, t, r)
494    elif diagn == 'LMDZrh':
495           
496        var0 = ncobj.variables[depvars[0]][:]
497        var1 = ncobj.variables[depvars[1]][:]
498        var2 = ncobj.variables[depvars[2]][:]
499
500        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnames,dvnames)
501        ncvar.insert_variable(ncobj, 'hur', diagout, diagoutd, diagoutvd, newnc)
502
503# LMDZrhs (psol, t2m, q2m)
504    elif diagn == 'LMDZrhs':
505           
506        var0 = ncobj.variables[depvars[0]][:]
507        var1 = ncobj.variables[depvars[1]][:]
508        var2 = ncobj.variables[depvars[2]][:]
509
510        dnamesvar = ncobj.variables[depvars[0]].dimensions
511        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
512
513        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
514
515        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
516
517# mrso: total soil moisture SMOIS, DZS
518    elif diagn == 'WRFmrso':
519           
520        var0 = ncobj.variables[depvars[0]][:]
521        var10 = ncobj.variables[depvars[1]][:]
522        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
523        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
524
525        var1 = var0.copy()*0.
526        var2 = var0.copy()*0.+1.
527        # Must be a better way....
528        for j in range(var0.shape[2]):
529          for i in range(var0.shape[3]):
530              var1[:,:,j,i] = var10
531
532        diagout, diagoutd, diagoutvd = diag.Forcompute_zint(var0, var1, var2,        \
533          dnamesvar, dvnamesvar)
534
535        # Removing the nonChecking variable-dimensions from the initial list
536        varsadd = []
537        diagoutvd = list(dvnames)
538        for nonvd in NONchkvardims:
539            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
540            varsadd.append(nonvd)
541        ncvar.insert_variable(ncobj, 'mrso', diagout, diagoutd, diagoutvd, newnc)
542
543# mslp: mean sea level pressure (pres, psfc, terrain, temp, qv)
544    elif diagn == 'mslp' or diagn == 'WRFmslp':
545           
546        var1 = ncobj.variables[depvars[1]][:]
547        var2 = ncobj.variables[depvars[2]][:]
548        var4 = ncobj.variables[depvars[4]][:]
549
550        if diagn == 'WRFmslp':
551            var0 = WRFp
552            var3 = WRFt
553            dnamesvar = ncobj.variables['P'].dimensions
554        else:
555            var0 = ncobj.variables[depvars[0]][:]
556            var3 = ncobj.variables[depvars[3]][:]
557            dnamesvar = ncobj.variables[depvars[0]].dimensions
558
559        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
560
561        diagout, diagoutd, diagoutvd = diag.compute_mslp(var0, var1, var2, var3, var4,    \
562          dnamesvar, dvnamesvar)
563
564        # Removing the nonChecking variable-dimensions from the initial list
565        varsadd = []
566        diagoutvd = list(dvnames)
567        for nonvd in NONchkvardims:
568            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
569            varsadd.append(nonvd)
570        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
571
572# OMEGAw (omega, p, t) from NCL formulation (https://www.ncl.ucar.edu/Document/Functions/Contributed/omega_to_w.shtml)
573    elif diagn == 'OMEGAw':
574           
575        var0 = ncobj.variables[depvars[0]][:]
576        var1 = ncobj.variables[depvars[1]][:]
577        var2 = ncobj.variables[depvars[2]][:]
578
579        dnamesvar = ncobj.variables[depvars[0]].dimensions
580        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
581
582        diagout, diagoutd, diagoutvd = diag.compute_OMEGAw(var0,var1,var2,dnamesvar,dvnamesvar)
583
584        ncvar.insert_variable(ncobj, 'wa', diagout, diagoutd, diagoutvd, newnc)
585
586# raintot: instantaneous total precipitation from WRF as (RAINC + RAINC) / dTime
587    elif diagn == 'RAINTOT':
588
589        var0 = ncobj.variables[depvars[0]]
590        var1 = ncobj.variables[depvars[1]]
591        if depvars[2] != 'WRFtime':
592            var2 = ncobj.variables[depvars[2]]
593        else:
594            var2 = np.arange(var0.shape[0], dtype=int)
595
596        var = var0[:] + var1[:]
597
598        dnamesvar = var0.dimensions
599        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
600
601        diagout, diagoutd, diagoutvd = diag.compute_deaccum(var,dnamesvar,dvnamesvar)
602
603# Transforming to a flux
604        if var2.shape[0] > 1:
605            if depvars[2] != 'WRFtime':
606                dtimeunits = var2.getncattr('units')
607                tunits = dtimeunits.split(' ')[0]
608   
609                dtime = (var2[1] - var2[0])*timeunits_seconds(tunits)
610            else:
611                var2 = ncobj.variables['Times']
612                time1 = var2[0,:]
613                time2 = var2[1,:]
614                tmf1 = ''
615                tmf2 = ''
616                for ic in range(len(time1)):
617                    tmf1 = tmf1 + time1[ic]
618                    tmf2 = tmf2 + time2[ic]
619                dtdate1 = dtime.datetime.strptime(tmf1,"%Y-%m-%d_%H:%M:%S")
620                dtdate2 = dtime.datetime.strptime(tmf2,"%Y-%m-%d_%H:%M:%S")
621                diffdate12 = dtdate2 - dtdate1
622                dtime = diffdate12.total_seconds()
623                print 'dtime:',dtime
624        else:
625            print warnmsg
626            print '  ' + main + ": only 1 time-step for '" + diag + "' !!"
627            print '    leaving a zero value!'
628            diagout = var0[:]*0.
629            dtime=1.
630
631        # Removing the nonChecking variable-dimensions from the initial list
632        varsadd = []
633        for nonvd in NONchkvardims:
634            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
635            varsadd.append(nonvd)
636           
637        ncvar.insert_variable(ncobj, 'pr', diagout/dtime, diagoutd, diagoutvd, newnc)
638
639# rhs (psfc, t, q) from TimeSeries files
640    elif diagn == 'TSrhs':
641           
642        p0=100000.
643        var0 = ncobj.variables[depvars[0]][:]
644        var1 = (ncobj.variables[depvars[1]][:])*(var0/p0)**(2./7.)
645        var2 = ncobj.variables[depvars[2]][:]
646
647        dnamesvar = ncobj.variables[depvars[0]].dimensions
648        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
649
650        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
651
652        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
653
654# slw: total soil liquid water SH2O, DZS
655    elif diagn == 'WRFslw':
656           
657        var0 = ncobj.variables[depvars[0]][:]
658        var10 = ncobj.variables[depvars[1]][:]
659        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
660        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
661
662        var1 = var0.copy()*0.
663        var2 = var0.copy()*0.+1.
664        # Must be a better way....
665        for j in range(var0.shape[2]):
666          for i in range(var0.shape[3]):
667              var1[:,:,j,i] = var10
668
669        diagout, diagoutd, diagoutvd = diag.Forcompute_zint(var0, var1, var2,        \
670          dnamesvar, dvnamesvar)
671
672        # Removing the nonChecking variable-dimensions from the initial list
673        varsadd = []
674        diagoutvd = list(dvnames)
675        for nonvd in NONchkvardims:
676            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
677            varsadd.append(nonvd)
678        ncvar.insert_variable(ncobj, 'slw', diagout, diagoutd, diagoutvd, newnc)
679
680# td (psfc, t, q) from TimeSeries files
681    elif diagn == 'TStd' or diagn == 'td':
682           
683        var0 = ncobj.variables[depvars[0]][:]
684        var1 = ncobj.variables[depvars[1]][:] - 273.15
685        var2 = ncobj.variables[depvars[2]][:]
686
687        dnamesvar = ncobj.variables[depvars[0]].dimensions
688        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
689
690        diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
691
692        ncvar.insert_variable(ncobj, 'tds', diagout, diagoutd, diagoutvd, newnc)
693
694# td (psfc, t, q) from TimeSeries files
695    elif diagn == 'TStdC' or diagn == 'tdC':
696           
697        var0 = ncobj.variables[depvars[0]][:]
698# Temperature is already in degrees Celsius
699        var1 = ncobj.variables[depvars[1]][:]
700        var2 = ncobj.variables[depvars[2]][:]
701
702        dnamesvar = ncobj.variables[depvars[0]].dimensions
703        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
704
705        diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
706
707        ncvar.insert_variable(ncobj, 'tds', diagout, diagoutd, diagoutvd, newnc)
708
709# wds (u, v)
710    elif diagn == 'TSwds' or diagn == 'wds' :
711 
712        var0 = ncobj.variables[depvars[0]][:]
713        var1 = ncobj.variables[depvars[1]][:]
714
715        dnamesvar = ncobj.variables[depvars[0]].dimensions
716        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
717
718        diagout, diagoutd, diagoutvd = diag.compute_wds(var0,var1,dnamesvar,dvnamesvar)
719
720        ncvar.insert_variable(ncobj, 'wds', diagout, diagoutd, diagoutvd, newnc)
721
722# wss (u, v)
723    elif diagn == 'TSwss' or diagn == 'wss':
724           
725        var0 = ncobj.variables[depvars[0]][:]
726        var1 = ncobj.variables[depvars[1]][:]
727
728        dnamesvar = ncobj.variables[depvars[0]].dimensions
729        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
730
731        diagout, diagoutd, diagoutvd = diag.compute_wss(var0,var1,dnamesvar,dvnamesvar)
732
733        ncvar.insert_variable(ncobj, 'wss', diagout, diagoutd, diagoutvd, newnc)
734
735# turbulence (var)
736    elif diagn == 'turbulence':
737
738        var0 = ncobj.variables[depvars][:]
739
740        dnamesvar = list(ncobj.variables[depvars].dimensions)
741        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
742
743        diagout, diagoutd, diagoutvd = diag.compute_turbulence(var0,dnamesvar,dvnamesvar)
744        valsvar = gen.variables_values(depvars)
745
746        newvarn = depvars + 'turb'
747        print main + '; Lluis newvarn:', newvarn
748        ncvar.insert_variable(ncobj, newvarn, diagout, diagoutd, 
749          diagoutvd, newnc)
750        print main + '; Lluis variables:', newnc.variables.keys()
751        varobj = newnc.variables[newvarn]
752        attrv = varobj.long_name
753        attr = varobj.delncattr('long_name')
754        newattr = ncvar.set_attribute(varobj, 'long_name', attrv +                   \
755          " Taylor decomposition turbulence term")
756
757# WRFbils fom WRF as HFX + LH
758    elif diagn == 'WRFbils':
759           
760        var0 = ncobj.variables[depvars[0]][:]
761        var1 = ncobj.variables[depvars[1]][:]
762
763        diagout = var0 + var1
764        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
765        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
766
767        ncvar.insert_variable(ncobj, 'bils', diagout, dnamesvar, dvnamesvar, newnc)
768
769# WRFcape_afwa CAPE, CIN, ZLFC, PLFC, LI following WRF 'phys/module_diaf_afwa.F'
770#   methodology as WRFt, WRFrh, WRFp, WRFgeop, HGT
771    elif diagn == 'WRFcape_afwa':
772        var0 = WRFt
773        var1 = WRFrh
774        var2 = WRFp
775        dz = WRFgeop.shape[1]
776        # de-staggering
777        var3 = 0.5*(WRFgeop[:,0:dz-1,:,:]+WRFgeop[:,1:dz,:,:])
778        var4 = ncobj.variables[depvars[4]][0,:,:]
779
780        dnamesvar = list(ncobj.variables['T'].dimensions)
781        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
782
783        diagout = np.zeros(var0.shape, dtype=np.float)
784        diagout1, diagout2, diagout3, diagout4, diagout5, diagoutd, diagoutvd =      \
785          diag.Forcompute_cape_afwa(var0, var1, var2, var3, var4, 3, dnamesvar,      \
786          dvnamesvar)
787
788        # Removing the nonChecking variable-dimensions from the initial list
789        varsadd = []
790        for nonvd in NONchkvardims:
791            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
792            varsadd.append(nonvd)
793
794        ncvar.insert_variable(ncobj, 'cape', diagout1, diagoutd, diagoutvd, newnc)
795        ncvar.insert_variable(ncobj, 'cin', diagout2, diagoutd, diagoutvd, newnc)
796        ncvar.insert_variable(ncobj, 'zlfc', diagout3, diagoutd, diagoutvd, newnc)
797        ncvar.insert_variable(ncobj, 'plfc', diagout4, diagoutd, diagoutvd, newnc)
798        ncvar.insert_variable(ncobj, 'li', diagout5, diagoutd, diagoutvd, newnc)
799
800# WRFclivi WRF water vapour path WRFdens, QICE, QGRAUPEL, QHAIL
801    elif diagn == 'WRFclivi':
802           
803        var0 = WRFdens
804        qtot = ncobj.variables[depvars[1]]
805        qtotv = qtot[:]
806        Nspecies = len(depvars) - 2
807        for iv in range(Nspecies):
808            if ncobj.variables.has_key(depvars[iv+2]):
809                var1 = ncobj.variables[depvars[iv+2]][:]
810                qtotv = qtotv + var1
811
812        dnamesvar = list(qtot.dimensions)
813        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
814
815        diagout, diagoutd, diagoutvd = diag.compute_clivi(var0, qtotv, dnamesvar,dvnamesvar)
816
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        ncvar.insert_variable(ncobj, 'clivi', diagout, diagoutd, diagoutvd, newnc)
824
825# WRFclwvi WRF water cloud-condensed path WRFdens, QCLOUD, QICE, QGRAUPEL, QHAIL
826    elif diagn == 'WRFclwvi':
827           
828        var0 = WRFdens
829        qtot = ncobj.variables[depvars[1]]
830        qtotv = ncobj.variables[depvars[1]]
831        Nspecies = len(depvars) - 2
832        for iv in range(Nspecies):
833            if ncobj.variables.has_key(depvars[iv+2]):
834                var1 = ncobj.variables[depvars[iv+2]]
835                qtotv = qtotv + var1[:]
836
837        dnamesvar = list(qtot.dimensions)
838        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
839
840        diagout, diagoutd, diagoutvd = diag.compute_clwvl(var0, qtotv, dnamesvar,dvnamesvar)
841
842        # Removing the nonChecking variable-dimensions from the initial list
843        varsadd = []
844        diagoutvd = list(dvnames)
845        for nonvd in NONchkvardims:
846            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
847            varsadd.append(nonvd)
848        ncvar.insert_variable(ncobj, 'clwvi', diagout, diagoutd, diagoutvd, newnc)
849
850# WRFgeop geopotential from WRF as PH + PHB
851    elif diagn == 'WRFgeop':
852        var0 = ncobj.variables[depvars[0]][:]
853        var1 = ncobj.variables[depvars[1]][:]
854
855        # de-staggering geopotential
856        diagout0 = var0 + var1
857        dt = diagout0.shape[0]
858        dz = diagout0.shape[1]
859        dy = diagout0.shape[2]
860        dx = diagout0.shape[3]
861
862        diagout = np.zeros((dt,dz-1,dy,dx), dtype=np.float)
863        diagout = 0.5*(diagout0[:,1:dz,:,:]+diagout0[:,0:dz-1,:,:])
864
865        # Removing the nonChecking variable-dimensions from the initial list
866        varsadd = []
867        diagoutvd = list(dvnames)
868        for nonvd in NONchkvardims:
869            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
870            varsadd.append(nonvd)
871
872        ncvar.insert_variable(ncobj, 'zg', diagout, dnames, diagoutvd, newnc)
873
874# WRFpotevap_orPM potential evapotranspiration following Penman-Monteith formulation
875#   implemented in ORCHIDEE as: WRFdens, UST, U10, V10, T2, PSFC, QVAPOR
876    elif diagn == 'WRFpotevap_orPM':
877        var0 = WRFdens[:,0,:,:]
878        var1 = ncobj.variables[depvars[1]][:]
879        var2 = ncobj.variables[depvars[2]][:]
880        var3 = ncobj.variables[depvars[3]][:]
881        var4 = ncobj.variables[depvars[4]][:]
882        var5 = ncobj.variables[depvars[5]][:]
883        var6 = ncobj.variables[depvars[6]][:,0,:,:]
884
885        dnamesvar = list(ncobj.variables[depvars[1]].dimensions)
886        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
887
888        diagout = np.zeros(var1.shape, dtype=np.float)
889        diagout, diagoutd, diagoutvd = diag.Forcompute_potevap_orPM(var0, var1, var2,\
890          var3, var4, var5, var6, dnamesvar, dvnamesvar)
891
892        # Removing the nonChecking variable-dimensions from the initial list
893        varsadd = []
894        for nonvd in NONchkvardims:
895            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
896            varsadd.append(nonvd)
897
898        ncvar.insert_variable(ncobj, 'evspsblpot', diagout, diagoutd, diagoutvd, newnc)
899
900# WRFmslp_ptarget sea-level pressure following ECMWF method as PSFC, HGT, WRFt, WRFp, ZNU, ZNW
901    elif diagn == 'WRFpsl_ecmwf':
902        var0 = ncobj.variables[depvars[0]][:]
903        var1 = ncobj.variables[depvars[1]][0,:,:]
904        var2 = WRFt[:,0,:,:]
905        var4 = WRFp[:,0,:,:]
906        var5 = ncobj.variables[depvars[4]][0,:]
907        var6 = ncobj.variables[depvars[5]][0,:]
908
909        # This is quite too appriximate!! passing pressure at half-levels to 2nd full
910        #   level, using eta values at full (ZNW) and half (ZNU) mass levels
911        var3 = WRFp[:,0,:,:] + (var6[1] - var5[0])*(WRFp[:,1,:,:] - WRFp[:,0,:,:])/  \
912          (var5[1]-var5[0])
913
914        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
915        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
916
917        diagout = np.zeros(var0.shape, dtype=np.float)
918        diagout, diagoutd, diagoutvd = diag.Forcompute_psl_ecmwf(var0, var1, var2,   \
919          var3, var4, dnamesvar, dvnamesvar)
920
921        # Removing the nonChecking variable-dimensions from the initial list
922        varsadd = []
923        for nonvd in NONchkvardims:
924            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
925            varsadd.append(nonvd)
926
927        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
928
929# WRFmslp_ptarget sea-level pressure following ptarget method as WRFp, PSFC, WRFt, HGT, QVAPOR
930    elif diagn == 'WRFpsl_ptarget':
931        var0 = WRFp
932        var1 = ncobj.variables[depvars[1]][:]
933        var2 = WRFt
934        var3 = ncobj.variables[depvars[3]][0,:,:]
935        var4 = ncobj.variables[depvars[4]][:]
936
937        dnamesvar = list(ncobj.variables[depvars[4]].dimensions)
938        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
939
940        diagout = np.zeros(var0.shape, dtype=np.float)
941        diagout, diagoutd, diagoutvd = diag.Forcompute_psl_ptarget(var0, var1, var2, \
942          var3, var4, 700000., dnamesvar, dvnamesvar)
943
944        # Removing the nonChecking variable-dimensions from the initial list
945        varsadd = []
946        for nonvd in NONchkvardims:
947            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
948            varsadd.append(nonvd)
949
950        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
951
952# WRFp pressure from WRF as P + PB
953    elif diagn == 'WRFp':
954           
955        diagout = WRFp
956
957        ncvar.insert_variable(ncobj, 'pres', diagout, dnames, dvnames, newnc)
958
959# WRFpos
960    elif diagn == 'WRFpos':
961           
962        dnamesvar = ncobj.variables['MAPFAC_M'].dimensions
963        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
964
965        ncvar.insert_variable(ncobj, 'WRFpos', WRFpos, dnamesvar, dvnamesvar, newnc)
966
967# WRFprw WRF water vapour path WRFdens, QVAPOR
968    elif diagn == 'WRFprw':
969           
970        var0 = WRFdens
971        var1 = ncobj.variables[depvars[1]]
972
973        dnamesvar = list(var1.dimensions)
974        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
975
976        diagout, diagoutd, diagoutvd = diag.compute_prw(var0, var1, dnamesvar,dvnamesvar)
977
978        # Removing the nonChecking variable-dimensions from the initial list
979        varsadd = []
980        diagoutvd = list(dvnames)
981        for nonvd in NONchkvardims:
982            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
983            varsadd.append(nonvd)
984        ncvar.insert_variable(ncobj, 'prw', diagout, diagoutd, diagoutvd, newnc)
985
986# WRFrh (P, T, QVAPOR)
987    elif diagn == 'WRFrh':
988           
989        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
990        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
991
992        ncvar.insert_variable(ncobj, 'hur', WRFrh, dnames, dvnames, newnc)
993
994# WRFrhs (PSFC, T2, Q2)
995    elif diagn == 'WRFrhs':
996           
997        var0 = ncobj.variables[depvars[0]][:]
998        var1 = ncobj.variables[depvars[1]][:]
999        var2 = ncobj.variables[depvars[2]][:]
1000
1001        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
1002        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1003
1004        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
1005        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
1006
1007# rvors (u10, v10, WRFpos)
1008    elif diagn == 'WRFrvors':
1009           
1010        var0 = ncobj.variables[depvars[0]]
1011        var1 = ncobj.variables[depvars[1]]
1012
1013        diagout = rotational_z(var0, var1, distx)
1014
1015        dnamesvar = ncobj.variables[depvars[0]].dimensions
1016        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1017
1018        ncvar.insert_variable(ncobj, 'rvors', diagout, dnamesvar, dvnamesvar, newnc)
1019
1020# WRFt (T, P, PB)
1021    elif diagn == 'WRFt':
1022        var0 = ncobj.variables[depvars[0]][:]
1023        var1 = ncobj.variables[depvars[1]][:]
1024        var2 = ncobj.variables[depvars[2]][:]
1025
1026        p0=100000.
1027        p=var1 + var2
1028
1029        WRFt = (var0 + 300.)*(p/p0)**(2./7.)
1030
1031        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1032        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1033
1034        # Removing the nonChecking variable-dimensions from the initial list
1035        varsadd = []
1036        diagoutvd = list(dvnames)
1037        for nonvd in NONchkvardims:
1038            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1039            varsadd.append(nonvd)
1040
1041        ncvar.insert_variable(ncobj, 'ta', WRFt, dnames, diagoutvd, newnc)
1042
1043# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
1044    elif diagn == 'WRFua':
1045        var0 = ncobj.variables[depvars[0]][:]
1046        var1 = ncobj.variables[depvars[1]][:]
1047        var2 = ncobj.variables[depvars[2]][:]
1048        var3 = ncobj.variables[depvars[3]][:]
1049
1050        # un-staggering variables
1051        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1052        ua = np.zeros(tuple(unstgdims), dtype=np.float)
1053        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1054        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1055        unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1056        unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1057
1058        for iz in range(var0.shape[1]):
1059            ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
1060
1061        dnamesvar = ['Time','bottom_top','south_north','west_east']
1062        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1063
1064        # Removing the nonChecking variable-dimensions from the initial list
1065        varsadd = []
1066        diagoutvd = list(dvnames)
1067        for nonvd in NONchkvardims:
1068            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1069            varsadd.append(nonvd)
1070
1071        ncvar.insert_variable(ncobj, 'ua', ua, dnames, diagoutvd, newnc)
1072
1073# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
1074    elif diagn == 'WRFva':
1075        var0 = ncobj.variables[depvars[0]][:]
1076        var1 = ncobj.variables[depvars[1]][:]
1077        var2 = ncobj.variables[depvars[2]][:]
1078        var3 = ncobj.variables[depvars[3]][:]
1079
1080        # un-staggering variables
1081        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1082        va = np.zeros(tuple(unstgdims), dtype=np.float)
1083        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1084        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1085        unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1086        unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1087        for iz in range(var0.shape[1]):
1088            va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
1089
1090        dnamesvar = ['Time','bottom_top','south_north','west_east']
1091        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1092
1093        # Removing the nonChecking variable-dimensions from the initial list
1094        varsadd = []
1095        diagoutvd = list(dvnames)
1096        for nonvd in NONchkvardims:
1097            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1098            varsadd.append(nonvd)
1099        ncvar.insert_variable(ncobj, 'va', va, dnames, diagoutvd, newnc)
1100
1101# WRFtime
1102    elif diagn == 'WRFtime':
1103           
1104        diagout = WRFtime
1105
1106        dnamesvar = ['Time']
1107        dvnamesvar = ['Times']
1108
1109        ncvar.insert_variable(ncobj, 'time', diagout, dnamesvar, dvnamesvar, newnc)
1110
1111# ws (U, V)
1112    elif diagn == 'ws':
1113           
1114        var0 = ncobj.variables[depvars[0]][:]
1115        var1 = ncobj.variables[depvars[1]][:]
1116        # un-staggering variables
1117        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1118        va = np.zeros(tuple(unstgdims), dtype=np.float)
1119        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1120        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1121        unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1122        unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1123
1124        dnamesvar = ['Time','bottom_top','south_north','west_east']
1125        diagout = np.sqrt(unstgvar0*unstgvar0 + unstgvar1*unstgvar1)
1126
1127#        dnamesvar = ncobj.variables[depvars[0]].dimensions
1128        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1129
1130        # Removing the nonChecking variable-dimensions from the initial list
1131        varsadd = []
1132        diagoutvd = list(dvnamesvar)
1133        for nonvd in NONchkvardims:
1134            if gen.searchInlist(dvnamesvar,nonvd): diagoutvd.remove(nonvd)
1135            varsadd.append(nonvd)
1136        ncvar.insert_variable(ncobj, 'ws', diagout, dnamesvar, diagoutvd, newnc)
1137
1138# wss (u10, v10)
1139    elif diagn == 'wss':
1140           
1141        var0 = ncobj.variables[depvars[0]][:]
1142        var1 = ncobj.variables[depvars[1]][:]
1143
1144        diagout = np.sqrt(var0*var0 + var1*var1)
1145
1146        dnamesvar = ncobj.variables[depvars[0]].dimensions
1147        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1148
1149        ncvar.insert_variable(ncobj, 'wss', diagout, dnamesvar, dvnamesvar, newnc)
1150
1151# WRFheight height from WRF geopotential as WRFGeop/g
1152    elif diagn == 'WRFheight':
1153           
1154        diagout = WRFgeop/grav
1155
1156        # Removing the nonChecking variable-dimensions from the initial list
1157        varsadd = []
1158        diagoutvd = list(dvnames)
1159        for nonvd in NONchkvardims:
1160            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1161            varsadd.append(nonvd)
1162
1163        ncvar.insert_variable(ncobj, 'zhgt', diagout, dnames, diagoutvd, newnc)
1164
1165# WRFheightrel relative-height from WRF geopotential as WRFgeop(PH + PHB)/g-HGT 'WRFheightrel|PH@PHB@HGT
1166    elif diagn == 'WRFheightrel':
1167        var0 = ncobj.variables[depvars[0]][:]
1168        var1 = ncobj.variables[depvars[1]][:]
1169        var2 = ncobj.variables[depvars[2]][:]
1170
1171        dimz = var0.shape[1]
1172        diagout = np.zeros(tuple(var0.shape), dtype=np.float)
1173        for iz in range(dimz):
1174            diagout[:,iz,:,:] = (var0[:,iz,:,:]+ var1[:,iz,:,:])/grav - var2
1175
1176        # Removing the nonChecking variable-dimensions from the initial list
1177        varsadd = []
1178        diagoutvd = list(dvnames)
1179        for nonvd in NONchkvardims:
1180            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1181            varsadd.append(nonvd)
1182
1183        ncvar.insert_variable(ncobj, 'zhgtrel', diagout, dnames, diagoutvd, newnc)
1184
1185# WRFzmla_gen generic boundary layer hieght computation from WRF theta, QVAPOR, WRFgeop, HGT,
1186    elif diagn == 'WRFzmlagen':
1187        var0 = ncobj.variables[depvars[0]][:]+300.
1188        var1 = ncobj.variables[depvars[1]][:]
1189        dimz = var0.shape[1]
1190        var2 = WRFgeop[:,1:dimz+1,:,:]/9.8
1191        var3 = ncobj.variables[depvars[3]][0,:,:]
1192
1193        diagout, diagoutd, diagoutvd = diag.Forcompute_zmla_gen(var0,var1,var2,var3, \
1194          dnames,dvnames)
1195
1196        # Removing the nonChecking variable-dimensions from the initial list
1197        varsadd = []
1198        for nonvd in NONchkvardims:
1199            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1200            varsadd.append(nonvd)
1201
1202        ncvar.insert_variable(ncobj, 'zmla', diagout, diagoutd, diagoutvd, newnc)
1203
1204# WRFzwind wind extrapolation at a given height using power law computation from WRF
1205#   U, V, WRFz, U10, V10, SINALPHA, COSALPHA, z=[zval]
1206    elif diagn == 'WRFzwind':
1207        var0 = ncobj.variables[depvars[0]][:]
1208        var1 = ncobj.variables[depvars[1]][:]
1209        var2 = WRFz
1210        var3 = ncobj.variables[depvars[3]][:]
1211        var4 = ncobj.variables[depvars[4]][:]
1212        var5 = ncobj.variables[depvars[5]][0,:,:]
1213        var6 = ncobj.variables[depvars[6]][0,:,:]
1214        var7 = np.float(depvars[7].split('=')[1])
1215
1216        # un-staggering 3D winds
1217        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1218        va = np.zeros(tuple(unstgdims), dtype=np.float)
1219        unvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1220        unvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1221        unvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1222        unvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1223
1224        diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwind(unvar0,      \
1225          unvar1, var2, var3, var4, var5, var6, var7, dnames, dvnames)
1226
1227        # Removing the nonChecking variable-dimensions from the initial list
1228        varsadd = []
1229        for nonvd in NONchkvardims:
1230            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1231            varsadd.append(nonvd)
1232
1233        ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
1234        ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)
1235
1236# WRFzwind wind extrapolation at a given hieght using logarithmic law computation
1237#   from WRF U, V, WRFz, U10, V10, SINALPHA, COSALPHA, z=[zval]
1238    elif diagn == 'WRFzwind_log':
1239        var0 = ncobj.variables[depvars[0]][:]
1240        var1 = ncobj.variables[depvars[1]][:]
1241        var2 = WRFz
1242        var3 = ncobj.variables[depvars[3]][:]
1243        var4 = ncobj.variables[depvars[4]][:]
1244        var5 = ncobj.variables[depvars[5]][0,:,:]
1245        var6 = ncobj.variables[depvars[6]][0,:,:]
1246        var7 = np.float(depvars[7].split('=')[1])
1247
1248        # un-staggering 3D winds
1249        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1250        va = np.zeros(tuple(unstgdims), dtype=np.float)
1251        unvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1252        unvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1253        unvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1254        unvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1255
1256        diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwind_log(unvar0,  \
1257          unvar1, var2, var3, var4, var5, var6, var7, dnames, dvnames)
1258
1259        # Removing the nonChecking variable-dimensions from the initial list
1260        varsadd = []
1261        for nonvd in NONchkvardims:
1262            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1263            varsadd.append(nonvd)
1264
1265        ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
1266        ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)
1267
1268# WRFzwindMO wind extrapolation at a given height computation using Monin-Obukhow
1269#   theory from WRF UST, ZNT, RMOL, U10, V10, SINALPHA, COSALPHA, z=[zval]
1270#   NOTE: only useful for [zval] < 80. m
1271    elif diagn == 'WRFzwindMO':
1272        var0 = ncobj.variables[depvars[0]][:]
1273        var1 = ncobj.variables[depvars[1]][:]
1274        var2 = ncobj.variables[depvars[2]][:]
1275        var3 = ncobj.variables[depvars[3]][:]
1276        var4 = ncobj.variables[depvars[4]][:]
1277        var5 = ncobj.variables[depvars[5]][0,:,:]
1278        var6 = ncobj.variables[depvars[6]][0,:,:]
1279        var7 = np.float(depvars[7].split('=')[1])
1280
1281        diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwindMO(var0, var1,\
1282          var2, var3, var4, var5, var6, var7, dnames, dvnames)
1283
1284        # Removing the nonChecking variable-dimensions from the initial list
1285        varsadd = []
1286        for nonvd in NONchkvardims:
1287            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1288            varsadd.append(nonvd)
1289
1290        ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
1291        ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)
1292
1293    else:
1294        print errormsg
1295        print '  ' + main + ": diagnostic '" + diagn + "' not ready!!!"
1296        print '    available diagnostics: ', availdiags
1297        quit(-1)
1298
1299    newnc.sync()
1300    # Adding that additional variables required to compute some diagnostics which
1301    #   where not in the original file
1302    for vadd in varsadd:
1303        if not gen.searchInlist(newnc.variables.keys(),vadd):
1304            attrs = dictcompvars[vadd]
1305            vvn = attrs['name']
1306            if not gen.searchInlist(newnc.variables.keys(), vvn):
1307                iidvn = dvnames.index(vadd)
1308                dnn = dnames[iidvn]
1309                if vadd == 'WRFtime':
1310                    dvarvals = WRFtime[:]
1311                newvar = newnc.createVariable(vvn, 'f8', (dnn))
1312                newvar[:] = dvarvals
1313                for attn in attrs.keys():
1314                    if attn != 'name':
1315                        attv = attrs[attn]
1316                        ncvar.set_attribute(newvar, attn, attv)
1317
1318#   end of diagnostics
1319
1320# Global attributes
1321##
1322ncvar.add_global_PyNCplot(newnc, main, None, '2.0')
1323
1324gorigattrs = ncobj.ncattrs()
1325for attr in gorigattrs:
1326    attrv = ncobj.getncattr(attr)
1327    atvar = ncvar.set_attribute(newnc, attr, attrv)
1328
1329ncobj.close()
1330newnc.close()
1331
1332print '\n' + main + ': successfull writting of diagnostics file "' + ofile + '" !!!'
Note: See TracBrowser for help on using the repository browser.