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

Last change on this file since 1624 was 1586, checked in by lfita, 8 years ago

Adding the removing of the non-checking variables for 'prw'
File size: 63.8 KB
RevLine 
[365]1# Pthon 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
[772]4#   To be used with module_ForDiagnostics.F90, module_ForDiagnosticsVars.F90, module_generic.F90
[1150]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
[1149]7
[413]8## e.g. # diagnostics.py -d 'Time@time,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
[365]10
11from optparse import OptionParser
12import numpy as np
13from netCDF4 import Dataset as NetCDFFile
14import os
15import re
16import nc_var_tools as ncvar
[756]17import generic_tools as gen
[654]18import datetime as dtime
[1163]19import module_ForDiag as fdin
[365]20
21main = 'diagnostics.py'
22errormsg = 'ERROR -- error -- ERROR -- error'
23warnmsg = 'WARNING -- warning -- WARNING -- warning'
24
[654]25# Constants
26grav = 9.81
27
[365]28# Gneral information
29##
30def reduce_spaces(string):
31    """ Function to give words of a line of text removing any extra space
32    """
33    values = string.replace('\n','').split(' ')
34    vals = []
35    for val in values:
36         if len(val) > 0:
37             vals.append(val)
38
39    return vals
40
41def variable_combo(varn,combofile):
42    """ Function to provide variables combination from a given variable name
43      varn= name of the variable
44      combofile= ASCII file with the combination of variables
45        [varn] [combo]
46          [combo]: '@' separated list of variables to use to generate [varn]
47            [WRFdt] to get WRF time-step (from general attributes)
48    >>> variable_combo('WRFprls','/home/lluis/PY/diagnostics.inf')
49    deaccum@RAINNC@XTIME@prnc
50    """
51    fname = 'variable_combo'
52
53    if varn == 'h':
54        print fname + '_____________________________________________________________'
55        print variable_combo.__doc__
56        quit()
57
58    if not os.path.isfile(combofile):
59        print errormsg
60        print '  ' + fname + ": file with combinations '" + combofile +              \
61          "' does not exist!!"
62        quit(-1)
63
64    objf = open(combofile, 'r')
65
66    found = False
67    for line in objf:
68        linevals = reduce_spaces(line)
69        varnf = linevals[0]
70        combo = linevals[1].replace('\n','')
71        if varn == varnf: 
72            found = True
73            break
74
75    if not found:
76        print errormsg
77        print '  ' + fname + ": variable '" + varn + "' not found in '" + combofile +\
78          "' !!"
79        combo='ERROR'
80
81    objf.close()
82
83    return combo
84
85# Mathematical operators
86##
[649]87def compute_accum(varv, dimns, dimvns):
88    """ Function to compute the accumulation of a variable
89    compute_accum(varv, dimnames, dimvns)
90      [varv]= values to accum (assuming [t,])
91      [dimns]= list of the name of the dimensions of the [varv]
92      [dimvns]= list of the name of the variables with the values of the
93        dimensions of [varv]
94    """
95    fname = 'compute_accum'
96
97    deacdims = dimns[:]
98    deacvdims = dimvns[:]
99
100    slicei = []
101    slicee = []
102
103    Ndims = len(varv.shape)
104    for iid in range(0,Ndims):
105        slicei.append(slice(0,varv.shape[iid]))
106        slicee.append(slice(0,varv.shape[iid]))
107
108    slicee[0] = np.arange(varv.shape[0])
109    slicei[0] = np.arange(varv.shape[0])
110    slicei[0][1:varv.shape[0]] = np.arange(varv.shape[0]-1)
111
112    vari = varv[tuple(slicei)]
113    vare = varv[tuple(slicee)]
114
115    ac = vari*0.
116    for it in range(1,varv.shape[0]):
117        ac[it,] = ac[it-1,] + vare[it,]
118
119    return ac, deacdims, deacvdims
120
[365]121def compute_deaccum(varv, dimns, dimvns):
122    """ Function to compute the deaccumulation of a variable
123    compute_deaccum(varv, dimnames, dimvns)
124      [varv]= values to deaccum (assuming [t,])
125      [dimns]= list of the name of the dimensions of the [varv]
126      [dimvns]= list of the name of the variables with the values of the
127        dimensions of [varv]
128    """
129    fname = 'compute_deaccum'
130
131    deacdims = dimns[:]
132    deacvdims = dimvns[:]
133
134    slicei = []
135    slicee = []
136
137    Ndims = len(varv.shape)
138    for iid in range(0,Ndims):
139        slicei.append(slice(0,varv.shape[iid]))
140        slicee.append(slice(0,varv.shape[iid]))
141
142    slicee[0] = np.arange(varv.shape[0])
143    slicei[0] = np.arange(varv.shape[0])
144    slicei[0][1:varv.shape[0]] = np.arange(varv.shape[0]-1)
145
146    vari = varv[tuple(slicei)]
147    vare = varv[tuple(slicee)]
148
149    deac = vare - vari
150
151    return deac, deacdims, deacvdims
152
153def derivate_centered(var,dim,dimv):
154    """ Function to compute the centered derivate of a given field
155      centered derivate(n) = (var(n-1) + var(n+1))/(2*dn).
156    [var]= variable
157    [dim]= which dimension to compute the derivate
158    [dimv]= dimension values (can be of different dimension of [var])
159    >>> derivate_centered(np.arange(16).reshape(4,4)*1.,1,1.)
160    [[  0.   1.   2.   0.]
161     [  0.   5.   6.   0.]
162     [  0.   9.  10.   0.]
163     [  0.  13.  14.   0.]]
164    """
165
166    fname = 'derivate_centered'
167   
168    vark = var.dtype
169
170    if hasattr(dimv, "__len__"):
171# Assuming that the last dimensions of var [..., N, M] are the same of dimv [N, M]
172        if len(var.shape) != len(dimv.shape):
173            dimvals = np.zeros((var.shape), dtype=vark)
174            if len(var.shape) - len(dimv.shape) == 1:
175                for iz in range(var.shape[0]):
176                    dimvals[iz,] = dimv
177            elif len(var.shape) - len(dimv.shape) == 2:
178                for it in range(var.shape[0]):
179                    for iz in range(var.shape[1]):
180                        dimvals[it,iz,] = dimv
181            else:
182                print errormsg
183                print '  ' + fname + ': dimension difference between variable',      \
184                  var.shape,'and variable with dimension values',dimv.shape,         \
185                  ' not ready !!!'
186                quit(-1)
187        else:
188            dimvals = dimv
189    else:
190# dimension values are identical everywhere!
191# from: http://stackoverflow.com/questions/16807011/python-how-to-identify-if-a-variable-is-an-array-or-a-scalar   
192        dimvals = np.ones((var.shape), dtype=vark)*dimv
193
194    derivate = np.zeros((var.shape), dtype=vark)
195    if dim > len(var.shape) - 1:
196        print errormsg
197        print '  ' + fname + ': dimension',dim,' too big for given variable of ' +   \
198          'shape:', var.shape,'!!!'
199        quit(-1)
200
201    slicebef = []
202    sliceaft = []
203    sliceder = []
204
205    for id in range(len(var.shape)):
206        if id == dim:
207            slicebef.append(slice(0,var.shape[id]-2))
208            sliceaft.append(slice(2,var.shape[id]))
209            sliceder.append(slice(1,var.shape[id]-1))
210        else:
211            slicebef.append(slice(0,var.shape[id]))
212            sliceaft.append(slice(0,var.shape[id]))
213            sliceder.append(slice(0,var.shape[id]))
214
215    if hasattr(dimv, "__len__"):
216        derivate[tuple(sliceder)] = (var[tuple(slicebef)] + var[tuple(sliceaft)])/   \
217          ((dimvals[tuple(sliceaft)] - dimvals[tuple(slicebef)]))
218        print (dimvals[tuple(sliceaft)] - dimvals[tuple(slicebef)])
219    else:
220        derivate[tuple(sliceder)] = (var[tuple(slicebef)] + var[tuple(sliceaft)])/   \
221          (2.*dimv)
222
223#    print 'before________'
224#    print var[tuple(slicebef)]
225
226#    print 'after________'
227#    print var[tuple(sliceaft)]
228
229    return derivate
230
231def rotational_z(Vx,Vy,pos):
232    """ z-component of the rotatinoal of horizontal vectorial field
233    \/ x (Vx,Vy,Vz) = \/xVy - \/yVx
234    [Vx]= Variable component x
235    [Vy]=  Variable component y
236    [pos]= poisition of the grid points
237    >>> rotational_z(np.arange(16).reshape(4,4)*1., np.arange(16).reshape(4,4)*1., 1.)
238    [[  0.   1.   2.   0.]
239     [ -4.   0.   0.  -7.]
240     [ -8.   0.   0. -11.]
241     [  0.  13.  14.   0.]]
242    """
243
244    fname =  'rotational_z'
245
246    ndims = len(Vx.shape)
247    rot1 = derivate_centered(Vy,ndims-1,pos)
248    rot2 = derivate_centered(Vx,ndims-2,pos)
249
250    rot = rot1 - rot2
251
252    return rot
253
254# Diagnostics
255##
256
257def var_clt(cfra):
258    """ Function to compute the total cloud fraction following 'newmicro.F90' from
259      LMDZ using 1D vertical column values
260      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
261    """
262    ZEPSEC=1.0E-12
263
264    fname = 'var_clt'
265
266    zclear = 1.
267    zcloud = 0.
268
269    dz = cfra.shape[0]
270    for iz in range(dz):
271        zclear =zclear*(1.-np.max([cfra[iz],zcloud]))/(1.-np.min([zcloud,1.-ZEPSEC]))
272        clt = 1. - zclear
273        zcloud = cfra[iz]
274
275    return clt
276
277def compute_clt(cldfra, dimns, dimvns):
278    """ Function to compute the total cloud fraction following 'newmicro.F90' from
279      LMDZ
280    compute_clt(cldfra, dimnames)
281      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
282      [dimns]= list of the name of the dimensions of [cldfra]
283      [dimvns]= list of the name of the variables with the values of the
284        dimensions of [cldfra]
285    """
286    fname = 'compute_clt'
287
288    cltdims = dimns[:]
289    cltvdims = dimvns[:]
290
291    if len(cldfra.shape) == 4:
292        clt = np.zeros((cldfra.shape[0],cldfra.shape[2],cldfra.shape[3]),            \
293          dtype=np.float)
294        dx = cldfra.shape[3]
295        dy = cldfra.shape[2]
296        dz = cldfra.shape[1]
297        dt = cldfra.shape[0]
298        cltdims.pop(1)
299        cltvdims.pop(1)
300
301        for it in range(dt):
302            for ix in range(dx):
303                for iy in range(dy):
304                    zclear = 1.
305                    zcloud = 0.
[1581]306                    gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
[365]307                    clt[it,iy,ix] = var_clt(cldfra[it,:,iy,ix])
308
309    else:
310        clt = np.zeros((cldfra.shape[1],cldfra.shape[2]), dtype=np.float)
311        dx = cldfra.shape[2]
312        dy = cldfra.shape[1]
313        dy = cldfra.shape[0]
314        cltdims.pop(0)
315        cltvdims.pop(0)
316        for ix in range(dx):
317            for iy in range(dy):
318                zclear = 1.
319                zcloud = 0.
[1581]320                gen.percendone(ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
[365]321                clt[iy,ix] = var_clt(cldfra[:,iy,ix])
322
323    return clt, cltdims, cltvdims
324
[772]325def Forcompute_clt(cldfra, dimns, dimvns):
326    """ Function to compute the total cloud fraction following 'newmicro.F90' from
327      LMDZ via a Fortran module
328    compute_clt(cldfra, dimnames)
329      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
330      [dimns]= list of the name of the dimensions of [cldfra]
331      [dimvns]= list of the name of the variables with the values of the
332        dimensions of [cldfra]
333    """
334    fname = 'Forcompute_clt'
335
336    cltdims = dimns[:]
337    cltvdims = dimvns[:]
338
339
340    if len(cldfra.shape) == 4:
341        clt = np.zeros((cldfra.shape[0],cldfra.shape[2],cldfra.shape[3]),            \
342          dtype=np.float)
343        dx = cldfra.shape[3]
344        dy = cldfra.shape[2]
345        dz = cldfra.shape[1]
346        dt = cldfra.shape[0]
347        cltdims.pop(1)
348        cltvdims.pop(1)
349
350        clt = fdin.module_fordiagnostics.compute_clt4d2(cldfra[:])
351
352    else:
353        clt = np.zeros((cldfra.shape[1],cldfra.shape[2]), dtype=np.float)
354        dx = cldfra.shape[2]
355        dy = cldfra.shape[1]
356        dy = cldfra.shape[0]
357        cltdims.pop(0)
358        cltvdims.pop(0)
359
360        clt = fdin.module_fordiagnostics.compute_clt3d1(cldfra[:])
361
362    return clt, cltdims, cltvdims
363
[365]364def var_cllmh(cfra, p):
365    """ Fcuntion to compute cllmh on a 1D column
366    """
367
368    fname = 'var_cllmh'
369
370    ZEPSEC =1.0E-12
371    prmhc = 440.*100.
372    prmlc = 680.*100.
373
374    zclearl = 1.
375    zcloudl = 0.
376    zclearm = 1.
377    zcloudm = 0.
378    zclearh = 1.
379    zcloudh = 0.
380
381    dvz = cfra.shape[0]
382
383    cllmh = np.ones((3), dtype=np.float)
384
385    for iz in range(dvz):
386        if p[iz] < prmhc:
387            cllmh[2] = cllmh[2]*(1.-np.max([cfra[iz], zcloudh]))/(1.-                \
388              np.min([zcloudh,1.-ZEPSEC]))
389            zcloudh = cfra[iz]
390        elif p[iz] >= prmhc and p[iz] < prmlc:
391            cllmh[1] = cllmh[1]*(1.-np.max([cfra[iz], zcloudm]))/(1.-                \
392              np.min([zcloudm,1.-ZEPSEC]))
393            zcloudm = cfra[iz]
394        elif p[iz] >= prmlc:
395            cllmh[0] = cllmh[0]*(1.-np.max([cfra[iz], zcloudl]))/(1.-                \
396              np.min([zcloudl,1.-ZEPSEC]))
397            zcloudl = cfra[iz]
398
399    cllmh = 1.- cllmh
400
401    return cllmh
402
[772]403def Forcompute_cllmh(cldfra, pres, dimns, dimvns):
404    """ Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ via Fortran subroutine
405    compute_clt(cldfra, pres, dimns, dimvns)
406      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
407      [pres] = pressure field
408      [dimns]= list of the name of the dimensions of [cldfra]
409      [dimvns]= list of the name of the variables with the values of the
410        dimensions of [cldfra]
411    """
412    fname = 'Forcompute_cllmh'
413
414    cllmhdims = dimns[:]
415    cllmhvdims = dimvns[:]
416
417    if len(cldfra.shape) == 4:
418        dx = cldfra.shape[3]
419        dy = cldfra.shape[2]
420        dz = cldfra.shape[1]
421        dt = cldfra.shape[0]
422        cllmhdims.pop(1)
423        cllmhvdims.pop(1)
424
425        cllmh = fdin.module_fordiagnostics.compute_cllmh4d2(cldfra[:], pres[:])
426       
427    else:
428        dx = cldfra.shape[2]
429        dy = cldfra.shape[1]
430        dz = cldfra.shape[0]
431        cllmhdims.pop(0)
432        cllmhvdims.pop(0)
433
434        cllmh = fdin.module_fordiagnostics.compute_cllmh3d1(cldfra[:], pres[:])
435
436    return cllmh, cllmhdims, cllmhvdims
437
[365]438def compute_cllmh(cldfra, pres, dimns, dimvns):
439    """ Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ
440    compute_clt(cldfra, pres, dimns, dimvns)
441      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
442      [pres] = pressure field
443      [dimns]= list of the name of the dimensions of [cldfra]
444      [dimvns]= list of the name of the variables with the values of the
445        dimensions of [cldfra]
446    """
447    fname = 'compute_cllmh'
448
449    cllmhdims = dimns[:]
450    cllmhvdims = dimvns[:]
451
452    if len(cldfra.shape) == 4:
453        dx = cldfra.shape[3]
454        dy = cldfra.shape[2]
455        dz = cldfra.shape[1]
456        dt = cldfra.shape[0]
457        cllmhdims.pop(1)
458        cllmhvdims.pop(1)
459
460        cllmh = np.ones(tuple([3, dt, dy, dx]), dtype=np.float)
461
462        for it in range(dt):
463            for ix in range(dx):
464                for iy in range(dy):
[1581]465                    gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
[365]466                    cllmh[:,it,iy,ix] = var_cllmh(cldfra[it,:,iy,ix], pres[it,:,iy,ix])
467       
468    else:
469        dx = cldfra.shape[2]
470        dy = cldfra.shape[1]
471        dz = cldfra.shape[0]
472        cllmhdims.pop(0)
473        cllmhvdims.pop(0)
474
475        cllmh = np.ones(tuple([3, dy, dx]), dtype=np.float)
476
477        for ix in range(dx):
478            for iy in range(dy):
[1581]479                gen.percendone(ix*dy + iy,dx*dy, 5, 'diagnosted')
[365]480                cllmh[:,iy,ix] = var_cllmh(cldfra[:,iy,ix], pres[:,iy,ix])
481
482    return cllmh, cllmhdims, cllmhvdims
483
[1581]484def compute_clivi(dens, qtot, dimns, dimvns):
485    """ Function to compute cloud-ice water path (clivi)
486      [dens] = density [in kgkg-1] (assuming [t],z,y,x)
487      [qtot] = added mixing ratio of all cloud-ice species in [kgkg-1] (assuming [t],z,y,x)
488      [dimns]= list of the name of the dimensions of [q]
489      [dimvns]= list of the name of the variables with the values of the
490        dimensions of [q]
491    """
492    fname = 'compute_clivi'
493
494    clividims = dimns[:]
495    clivivdims = dimvns[:]
496
497    if len(qtot.shape) == 4:
498        clividims.pop(1)
499        clivivdims.pop(1)
500    else:
501        clividims.pop(0)
502        clivivdims.pop(0)
503
504    data1 = dens*qtot
505    clivi = np.sum(data1, axis=1)
506
507    return clivi, clividims, clivivdims
508
509
510def compute_clwvl(dens, qtot, dimns, dimvns):
511    """ Function to compute condensed water path (clwvl)
512      [dens] = density [in kgkg-1] (assuming [t],z,y,x)
513      [qtot] = added mixing ratio of all cloud-water species in [kgkg-1] (assuming [t],z,y,x)
514      [dimns]= list of the name of the dimensions of [q]
515      [dimvns]= list of the name of the variables with the values of the
516        dimensions of [q]
517    """
518    fname = 'compute_clwvl'
519
520    clwvldims = dimns[:]
521    clwvlvdims = dimvns[:]
522
523    if len(qtot.shape) == 4:
524        clwvldims.pop(1)
525        clwvlvdims.pop(1)
526    else:
527        clwvldims.pop(0)
528        clwvlvdims.pop(0)
529
530    data1 = dens*qtot
531    clwvl = np.sum(data1, axis=1)
532
533    return clwvl, clwvldims, clwvlvdims
534
[365]535def var_virtualTemp (temp,rmix):
536    """ This function returns virtual temperature in K,
537      temp: temperature [K]
538      rmix: mixing ratio in [kgkg-1]
539    """
540
541    fname = 'var_virtualTemp'
542
543    virtual=temp*(0.622+rmix)/(0.622*(1.+rmix))
544
545    return virtual
546
547
548def var_mslp(pres, psfc, ter, tk, qv):
549    """ Function to compute mslp on a 1D column
550    """
551
552    fname = 'var_mslp'
553
554    N = 1.0
555    expon=287.04*.0065/9.81
556    pref = 40000.
557
558# First find where about 400 hPa is located
559    dz=len(pres) 
560
561    kref = -1
562    pinc = pres[0] - pres[dz-1]
563
564    if pinc < 0.:
565        for iz in range(1,dz):
566            if pres[iz-1] >= pref and pres[iz] < pref: 
567                kref = iz
568                break
569    else:
570        for iz in range(dz-1):
571            if pres[iz] >= pref and pres[iz+1] < pref: 
572                kref = iz
573                break
574
575    if kref == -1:
576        print errormsg
577        print '  ' + fname + ': no reference pressure:',pref,'found!!'
578        print '    values:',pres[:]
579        quit(-1)
580
581    mslp = 0.
582
583# We are below both the ground and the lowest data level.
584
585# First, find the model level that is closest to a "target" pressure
586# level, where the "target" pressure is delta-p less that the local
587# value of a horizontally smoothed surface pressure field.  We use
588# delta-p = 150 hPa here. A standard lapse rate temperature profile
589# passing through the temperature at this model level will be used
590# to define the temperature profile below ground.  This is similar
591# to the Benjamin and Miller (1990) method, using 
592# 700 hPa everywhere for the "target" pressure.
593
594# ptarget = psfc - 15000.
595    ptarget = 70000.
596    dpmin=1.e4
597    kupper = 0
598    if pinc > 0.:
599        for iz in range(dz-1,0,-1):
600            kupper = iz
601            dp=np.abs( pres[iz] - ptarget )
602            if dp < dpmin: exit
603            dpmin = np.min([dpmin, dp])
604    else:
605        for iz in range(dz):
606            kupper = iz
607            dp=np.abs( pres[iz] - ptarget )
608            if dp < dpmin: exit
609            dpmin = np.min([dpmin, dp])
610
611    pbot=np.max([pres[0], psfc])
612#    zbot=0.
613
614#    tbotextrap=tk(i,j,kupper,itt)*(pbot/pres_field(i,j,kupper,itt))**expon
615#    tvbotextrap=virtual(tbotextrap,qv(i,j,1,itt))
616
617#    data_out(i,j,itt,1) = (zbot+tvbotextrap/.0065*(1.-(interp_levels(1)/pbot)**expon))
618    tbotextrap = tk[kupper]*(psfc/ptarget)**expon
619    tvbotextrap = var_virtualTemp(tbotextrap, qv[kupper])
620    mslp = psfc*( (tvbotextrap+0.0065*ter)/tvbotextrap)**(1./expon)
621
622    return mslp
623
624def compute_mslp(pressure, psurface, terrain, temperature, qvapor, dimns, dimvns):
625    """ Function to compute mslp: mean sea level pressure following p_interp.F90 from WRF
626    var_mslp(pres, ter, tk, qv, dimns, dimvns)
627      [pressure]= pressure field [Pa] (assuming [[t],z,y,x])
628      [psurface]= surface pressure field [Pa]
629      [terrain]= topography [m]
630      [temperature]= temperature [K]
631      [qvapor]= water vapour mixing ratio [kgkg-1]
632      [dimns]= list of the name of the dimensions of [cldfra]
633      [dimvns]= list of the name of the variables with the values of the
634        dimensions of [pres]
635    """
636
637    fname = 'compute_mslp'
638
639    mslpdims = list(dimns[:])
640    mslpvdims = list(dimvns[:])
641
642    if len(pressure.shape) == 4:
643        mslpdims.pop(1)
644        mslpvdims.pop(1)
645    else:
646        mslpdims.pop(0)
647        mslpvdims.pop(0)
648
649    if len(pressure.shape) == 4:
650        dx = pressure.shape[3]
651        dy = pressure.shape[2]
652        dz = pressure.shape[1]
653        dt = pressure.shape[0]
654
655        mslpv = np.zeros(tuple([dt, dy, dx]), dtype=np.float)
656
657# Terrain... to 2D !
658        terval = np.zeros(tuple([dy, dx]), dtype=np.float)
659        if len(terrain.shape) == 3:
660            terval = terrain[0,:,:]
661        else:
662            terval = terrain
663
664        for ix in range(dx):
665            for iy in range(dy):
666                if terval[iy,ix] > 0.:
667                    for it in range(dt):
668                        mslpv[it,iy,ix] = var_mslp(pressure[it,:,iy,ix],             \
669                          psurface[it,iy,ix], terval[iy,ix], temperature[it,:,iy,ix],\
670                          qvapor[it,:,iy,ix])
671
[1581]672                        gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
[365]673                else:
674                    mslpv[:,iy,ix] = psurface[:,iy,ix]
675
676    else:
677        dx = pressure.shape[2]
678        dy = pressure.shape[1]
679        dz = pressure.shape[0]
680
681        mslpv = np.zeros(tuple([dy, dx]), dtype=np.float)
682
683# Terrain... to 2D !
684        terval = np.zeros(tuple([dy, dx]), dtype=np.float)
685        if len(terrain.shape) == 3:
686            terval = terrain[0,:,:]
687        else:
688            terval = terrain
689
690        for ix in range(dx):
691            for iy in range(dy):
[1581]692                gen.percendone(ix*dy + iy,dx*dy, 5, 'diagnosted')
[365]693                if terval[iy,ix] > 0.:
694                    mslpv[iy,ix] = var_mslp(pressure[:,iy,ix], psurface[iy,ix],          \
695                      terval[iy,ix], temperature[:,iy,ix], qvapor[:,iy,ix])
696                else:
697                    mslpv[iy,ix] = psfc[iy,ix]
698
699    return mslpv, mslpdims, mslpvdims
700
[642]701def compute_OMEGAw(omega, p, t, dimns, dimvns):
702    """ Function to transform OMEGA [Pas-1] to velocities [ms-1]
703      tacking: https://www.ncl.ucar.edu/Document/Functions/Contributed/omega_to_w.shtml
704      [omega] = vertical velocity [in ms-1] (assuming [t],z,y,x)
705      [p] = pressure in [Pa] (assuming [t],z,y,x)
706      [t] = temperature in [K] (assuming [t],z,y,x)
707      [dimns]= list of the name of the dimensions of [q]
708      [dimvns]= list of the name of the variables with the values of the
709        dimensions of [q]
710    """
711    fname = 'compute_OMEGAw'
712
713    rgas = 287.058            # J/(kg-K) => m2/(s2 K)
714    g    = 9.80665            # m/s2
715
716    wdims = dimns[:]
717    wvdims = dimvns[:]
718
719    rho  = p/(rgas*t)         # density => kg/m3
720    w    = -omega/(rho*g)     
721
722    return w, wdims, wvdims
723
[365]724def compute_prw(dens, q, dimns, dimvns):
725    """ Function to compute water vapour path (prw)
726      [dens] = density [in kgkg-1] (assuming [t],z,y,x)
727      [q] = mixing ratio in [kgkg-1] (assuming [t],z,y,x)
728      [dimns]= list of the name of the dimensions of [q]
729      [dimvns]= list of the name of the variables with the values of the
730        dimensions of [q]
731    """
732    fname = 'compute_prw'
733
734    prwdims = dimns[:]
735    prwvdims = dimvns[:]
736
737    if len(q.shape) == 4:
738        prwdims.pop(1)
739        prwvdims.pop(1)
740    else:
741        prwdims.pop(0)
742        prwvdims.pop(0)
743
744    data1 = dens*q
745    prw = np.sum(data1, axis=1)
746
747    return prw, prwdims, prwvdims
748
749def compute_rh(p, t, q, dimns, dimvns):
750    """ Function to compute relative humidity following 'Tetens' equation (T,P) ...'
751      [t]= temperature (assuming [[t],z,y,x] in [K])
752      [p] = pressure field (assuming in [hPa])
753      [q] = mixing ratio in [kgkg-1]
754      [dimns]= list of the name of the dimensions of [t]
755      [dimvns]= list of the name of the variables with the values of the
756        dimensions of [t]
757    """
758    fname = 'compute_rh'
759
760    rhdims = dimns[:]
761    rhvdims = dimvns[:]
762
763    data1 = 10.*0.6112*np.exp(17.67*(t-273.16)/(t-29.65))
764    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
765
766    rh = q/data2
767
768    return rh, rhdims, rhvdims
769
[612]770def compute_td(p, temp, qv, dimns, dimvns):
771    """ Function to compute the dew point temperature
772      [p]= pressure [Pa]
773      [temp]= temperature [C]
774      [qv]= mixing ratio [kgkg-1]
775      [dimns]= list of the name of the dimensions of [p]
776      [dimvns]= list of the name of the variables with the values of the
777        dimensions of [p]
778    """
779    fname = 'compute_td'
780
781#    print '    ' + fname + ': computing dew-point temperature from TS as t and Tetens...'
782# tacking from: http://en.wikipedia.org/wiki/Dew_point
783    tk = temp
784    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
785    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
786
787    rh = qv/data2
788               
789    pa = rh * data1
[614]790    td = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
[612]791
792    tddims = dimns[:]
793    tdvdims = dimvns[:]
794
795    return td, tddims, tdvdims
796
[365]797def turbulence_var(varv, dimvn, dimn):
798    """ Function to compute the Taylor's decomposition turbulence term from a a given variable
799      x*=<x^2>_t-(<X>_t)^2
800    turbulence_var(varv,dimn)
801      varv= values of the variable
802      dimvn= names of the dimension of the variable
803      dimn= names of the dimensions (as a dictionary with 'X', 'Y', 'Z', 'T')
804    >>> turbulence_var(np.arange((27)).reshape(3,3,3),['time','y','x'],{'T':'time', 'Y':'y', 'X':'x'})
805    [[ 54.  54.  54.]
806     [ 54.  54.  54.]
807     [ 54.  54.  54.]]
808    """
809    fname = 'turbulence_varv'
810
811    timedimid = dimvn.index(dimn['T'])
812
813    varv2 = varv*varv
814
815    vartmean = np.mean(varv, axis=timedimid)
816    var2tmean = np.mean(varv2, axis=timedimid)
817
818    varvturb = var2tmean - (vartmean*vartmean)
819
820    return varvturb
821
822def compute_turbulence(v, dimns, dimvns):
823    """ Function to compute the rubulence term of the Taylor's decomposition ...'
824      x*=<x^2>_t-(<X>_t)^2
825      [v]= variable (assuming [[t],z,y,x])
826      [dimns]= list of the name of the dimensions of [v]
827      [dimvns]= list of the name of the variables with the values of the
828        dimensions of [v]
829    """
830    fname = 'compute_turbulence'
831
832    turbdims = dimns[:]
833    turbvdims = dimvns[:]
834
835    turbdims.pop(0)
836    turbvdims.pop(0)
837
838    v2 = v*v
839
840    vartmean = np.mean(v, axis=0)
841    var2tmean = np.mean(v2, axis=0)
842
843    turb = var2tmean - (vartmean*vartmean)
844
845    return turb, turbdims, turbvdims
846
[612]847def compute_wds(u, v, dimns, dimvns):
848    """ Function to compute the wind direction
849      [u]= W-E wind direction [ms-1, knot, ...]
850      [v]= N-S wind direction [ms-1, knot, ...]
851      [dimns]= list of the name of the dimensions of [u]
852      [dimvns]= list of the name of the variables with the values of the
853        dimensions of [u]
854    """
855    fname = 'compute_wds'
856
857#    print '    ' + fname + ': computing wind direction as ATAN2(v,u) ...'
858    theta = np.arctan2(v,u)
859    theta = np.where(theta < 0., theta + 2.*np.pi, theta)
860
861    wds = 360.*theta/(2.*np.pi)
862
863    wdsdims = dimns[:]
864    wdsvdims = dimvns[:]
865
866    return wds, wdsdims, wdsvdims
867
868def compute_wss(u, v, dimns, dimvns):
869    """ Function to compute the wind speed
870      [u]= W-E wind direction [ms-1, knot, ...]
871      [v]= N-S wind direction [ms-1, knot, ...]
872      [dimns]= list of the name of the dimensions of [u]
873      [dimvns]= list of the name of the variables with the values of the
874        dimensions of [u]
875    """
876    fname = 'compute_wss'
877
878#    print '    ' + fname + ': computing wind speed as SQRT(v**2 + u**2) ...'
879    wss = np.sqrt(u*u + v*v)
880
881    wssdims = dimns[:]
882    wssvdims = dimvns[:]
883
884    return wss, wssdims, wssvdims
885
[365]886def timeunits_seconds(dtu):
887    """ Function to transform a time units to seconds
888    timeunits_seconds(timeuv)
889      [dtu]= time units value to transform in seconds
890    """
891    fname='timunits_seconds'
892
893    if dtu == 'years':
894        times = 365.*24.*3600.
895    elif dtu == 'weeks':
896        times = 7.*24.*3600.
897    elif dtu == 'days':
898        times = 24.*3600.
899    elif dtu == 'hours':
900        times = 3600.
901    elif dtu == 'minutes':
902        times = 60.
903    elif dtu == 'seconds':
904        times = 1.
905    elif dtu == 'miliseconds':
906        times = 1./1000.
907    else:
908        print errormsg
909        print '  ' + fname  + ": time units '" + dtu + "' not ready !!"
910        quit(-1)
911
912    return times
913
914####### ###### ##### #### ### ## #
915comboinf="\nIF -d 'variable_combo', provides information of the combination to obtain -v [varn] with the ASCII file with the combinations as -f [combofile]"
916
917parser = OptionParser()
918parser.add_option("-f", "--netCDF_file", dest="ncfile", help="file to use", metavar="FILE")
919parser.add_option("-d", "--dimensions", dest="dimns", 
[1351]920  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)" + comboinf, 
[365]921  metavar="LABELS")
922parser.add_option("-v", "--variables", dest="varns", 
923  help=" [varn1]|[var11]@[...[varN1]],[...,[varnM]|[var1M]@[...[varLM]]] ',' list of variables to compute [varnK] and its necessary ones [var1K]...[varPK]", metavar="VALUES")
924
925(opts, args) = parser.parse_args()
926
927#######    #######
928## MAIN
929    #######
[649]930availdiags = ['ACRAINTOT', 'accum', 'clt', 'cllmh', 'deaccum', 'LMDZrh', 'mslp',     \
[1351]931  'OMEGAw', 'RAINTOT',                                                               \
[1581]932  'rvors', 'td', 'turbulence', 'WRFclivi', 'WRFclwvl', 'WRFgeop', 'WRFp',            \
933  'WRFrvors', 'ws', 'wds', 'wss', 'WRFheight', 'WRFheightrel', 'WRFua', 'WRFva']
[365]934
[649]935methods = ['accum', 'deaccum']
936
[365]937# Variables not to check
[612]938NONcheckingvars = ['cllmh', 'deaccum', 'TSrhs', 'TStd', 'TSwds', 'TSwss', 'WRFbils', \
[1581]939  'WRFclivi', 'WRFclwvl', 'WRFdens', 'WRFgeop',                                      \
[612]940  'WRFp', 'WRFtd',                                                                   \
[365]941  'WRFpos', 'WRFprc', 'WRFprls', 'WRFrh', 'LMDZrh', 'LMDZrhs', 'WRFrhs', 'WRFrvors', \
[914]942  'WRFt', 'WRFtime', 'WRFua', 'WRFva', 'WRFwds', 'WRFwss', 'WRFheight']
[365]943
[1351]944NONchkvardims = ['WRFtime']
945
[365]946ofile = 'diagnostics.nc'
947
948dimns = opts.dimns
949varns = opts.varns
950
951# Special method. knowing variable combination
952##
953if opts.dimns == 'variable_combo':
954    print warnmsg
955    print '  ' + main + ': knowing variable combination !!!'
956    combination = variable_combo(opts.varns,opts.ncfile)
957    print '     COMBO: ' + combination
958    quit(-1)
959
960if not os.path.isfile(opts.ncfile):
961    print errormsg
962    print '   ' + main + ": file '" + opts.ncfile + "' does not exist !!"
963    quit(-1)
964
965ncobj = NetCDFFile(opts.ncfile, 'r')
966
[1351]967# Looking for specific variables that might be use in more than one diagnostic
968WRFgeop_compute = False
969WRFp_compute = False
970WRFt_compute = False
971WRFrh_compute = False
972WRFght_compute = False
973WRFdens_compute = False
974WRFpos_compute = False
975WRFtime_compute = False
976
[365]977# File creation
978newnc = NetCDFFile(ofile,'w')
979
980# dimensions
981dimvalues = dimns.split(',')
982dnames = []
983dvnames = []
984
985for dimval in dimvalues:
[1351]986    dn = dimval.split('@')[0]
987    dnv = dimval.split('@')[1]
988    dnames.append(dn)
989    dvnames.append(dnv)
990    # Is there any dimension-variable which should be computed?
991    if dnv == 'WRFgeop':WRFgeop_compute = True
992    if dnv == 'WRFp': WRFp_compute = True
993    if dnv == 'WRFt': WRFt_compute = True
994    if dnv == 'WRFrh': WRFrh_compute = True
995    if dnv == 'WRFght': WRFght_compute = True
996    if dnv == 'WRFdens': WRFdens_compute = True
997    if dnv == 'WRFpos': WRFpos_compute = True
998    if dnv == 'WRFtime': WRFtime_compute = True
[365]999
1000# diagnostics to compute
1001diags = varns.split(',')
1002Ndiags = len(diags)
1003
1004for idiag in range(Ndiags):
1005    if diags[idiag].split('|')[1].find('@') == -1:
1006        depvars = diags[idiag].split('|')[1]
[654]1007        if depvars == 'WRFgeop':WRFgeop_compute = True
[365]1008        if depvars == 'WRFp': WRFp_compute = True
1009        if depvars == 'WRFt': WRFt_compute = True
1010        if depvars == 'WRFrh': WRFrh_compute = True
1011        if depvars == 'WRFght': WRFght_compute = True
1012        if depvars == 'WRFdens': WRFdens_compute = True
1013        if depvars == 'WRFpos': WRFpos_compute = True
[654]1014        if depvars == 'WRFtime': WRFtime_compute = True
[365]1015    else:
1016        depvars = diags[idiag].split('|')[1].split('@')
[756]1017        if gen.searchInlist(depvars, 'WRFgeop'): WRFgeop_compute = True
1018        if gen.searchInlist(depvars, 'WRFp'): WRFp_compute = True
1019        if gen.searchInlist(depvars, 'WRFt'): WRFt_compute = True
1020        if gen.searchInlist(depvars, 'WRFrh'): WRFrh_compute = True
1021        if gen.searchInlist(depvars, 'WRFght'): WRFght_compute = True
1022        if gen.searchInlist(depvars, 'WRFdens'): WRFdens_compute = True
1023        if gen.searchInlist(depvars, 'WRFpos'): WRFpos_compute = True
1024        if gen.searchInlist(depvars, 'WRFtime'): WRFtime_compute = True
[365]1025
[1351]1026# Dictionary with the new computed variables to be able to add them
1027dictcompvars = {}
[654]1028if WRFgeop_compute:
1029    print '  ' + main + ': Retrieving geopotential value from WRF as PH + PHB'
1030    dimv = ncobj.variables['PH'].shape
1031    WRFgeop = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
1032
[1351]1033    # Attributes of the variable
[1412]1034    Vvals = gen.variables_values('WRFgeop')
[1351]1035    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
1036      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
1037
[365]1038if WRFp_compute:
1039    print '  ' + main + ': Retrieving pressure value from WRF as P + PB'
1040    dimv = ncobj.variables['P'].shape
1041    WRFp = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
1042
[1351]1043    # Attributes of the variable
1044    Vvals = gen.variables_values('WRFp')
1045    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
1046      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
1047
[365]1048if WRFght_compute:
1049    print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...' 
1050    WRFght = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
1051
[1351]1052    # Attributes of the variable
1053    Vvals = gen.variables_values('WRFght')
1054    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
1055      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
1056
[365]1057if WRFrh_compute:
1058    print '    ' + main + ": computing relative humidity from WRF as 'Tetens'" +     \
1059      ' equation (T,P) ...'
1060    p0=100000.
1061    p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
1062    tk = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
1063    qv = ncobj.variables['QVAPOR'][:]
1064
1065    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
1066    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
1067
1068    WRFrh = qv/data2
1069
[1351]1070    # Attributes of the variable
1071    Vvals = gen.variables_values('WRFrh')
1072    dictcompvars['WRFrh'] = {'name': Vvals[0], 'standard_name': Vvals[1],            \
1073      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
1074
[365]1075if WRFt_compute:
1076    print '    ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
1077    p0=100000.
1078    p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
1079
1080    WRFt = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
1081
[1351]1082    # Attributes of the variable
1083    Vvals = gen.variables_values('WRFt')
1084    dictcompvars['WRFt'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
1085      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
1086
[365]1087if WRFdens_compute:
1088    print '    ' + main + ': computing air density from WRF as ((MU + MUB) * ' +     \
1089      'DNW)/g ...'
1090
1091# Just we need in in absolute values: Size of the central grid cell
1092##    dxval = ncobj.getncattr('DX')
1093##    dyval = ncobj.getncattr('DY')
1094##    mapfac = ncobj.variables['MAPFAC_M'][:]
1095##    area = dxval*dyval*mapfac
1096
1097    mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
1098    dnw = ncobj.variables['DNW'][:]
1099
1100    WRFdens = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]),        \
1101      dtype=np.float)
1102    levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)
1103
1104    for it in range(mu.shape[0]):
1105        for iz in range(dnw.shape[1]):
1106            levval.fill(np.abs(dnw[it,iz]))
1107            WRFdens[it,iz,:,:] = levval
1108            WRFdens[it,iz,:,:] = mu[it,:,:]*WRFdens[it,iz,:,:]/grav
1109
[1351]1110    # Attributes of the variable
1111    Vvals = gen.variables_values('WRFdens')
1112    dictcompvars['WRFdens'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
1113      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
1114
[365]1115if WRFpos_compute:
1116# WRF positions from the lowest-leftest corner of the matrix
1117    print '    ' + main + ': computing position from MAPFAC_M as sqrt(DY*j**2 + ' +  \
1118      'DX*x**2)*MAPFAC_M ...'
1119
1120    mapfac = ncobj.variables['MAPFAC_M'][:]
1121
1122    distx = np.float(ncobj.getncattr('DX'))
1123    disty = np.float(ncobj.getncattr('DY'))
1124
1125    print 'distx:',distx,'disty:',disty
1126
1127    dx = mapfac.shape[2]
1128    dy = mapfac.shape[1]
1129    dt = mapfac.shape[0]
1130
1131    WRFpos = np.zeros((dt, dy, dx), dtype=np.float)
1132
1133    for i in range(1,dx):
1134        WRFpos[0,0,i] = distx*i/mapfac[0,0,i]
1135    for j in range(1,dy):
1136        i=0
1137        WRFpos[0,j,i] = WRFpos[0,j-1,i] + disty/mapfac[0,j,i]
1138        for i in range(1,dx):
1139#            WRFpos[0,j,i] = np.sqrt((disty*j)**2. + (distx*i)**2.)/mapfac[0,j,i]
1140#            WRFpos[0,j,i] = np.sqrt((disty*j)**2. + (distx*i)**2.)
1141             WRFpos[0,j,i] = WRFpos[0,j,i-1] + distx/mapfac[0,j,i]
1142
1143    for it in range(1,dt):
1144        WRFpos[it,:,:] = WRFpos[0,:,:]
1145
[654]1146if WRFtime_compute:
1147    print '    ' + main + ': computing time from WRF as CFtime(Times) ...'
1148
1149    refdate='19491201000000'
1150    tunitsval='minutes'
1151
1152    timeobj = ncobj.variables['Times']
1153    timewrfv = timeobj[:]
1154
1155    yrref=refdate[0:4]
1156    monref=refdate[4:6]
1157    dayref=refdate[6:8]
1158    horref=refdate[8:10]
1159    minref=refdate[10:12]
1160    secref=refdate[12:14]
1161
1162    refdateS = yrref + '-' + monref + '-' + dayref + ' ' + horref + ':' + minref +   \
1163      ':' + secref
1164
1165    dt = timeobj.shape[0]
1166    WRFtime = np.zeros((dt), dtype=np.float)
1167
1168    for it in range(dt):
[865]1169        wrfdates = gen.datetimeStr_conversion(timewrfv[it,:],'WRFdatetime', 'matYmdHMS')
1170        WRFtime[it] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
[654]1171
1172    tunits = tunitsval + ' since ' + refdateS
1173
[1351]1174    # Attributes of the variable
1175    dictcompvars['WRFtime'] = {'name': 'time', 'standard_name': 'time',              \
1176      'long_name': 'time', 'units': tunits, 'calendar': 'gregorian'}
1177
[365]1178### ## #
1179# Going for the diagnostics
1180### ## #
1181print '  ' + main + ' ...'
[1404]1182varsadd = []
[365]1183
1184for idiag in range(Ndiags):
1185    print '    diagnostic:',diags[idiag]
1186    diag = diags[idiag].split('|')[0]
1187    depvars = diags[idiag].split('|')[1].split('@')
1188    if diags[idiag].split('|')[1].find('@') != -1:
1189        depvars = diags[idiag].split('|')[1].split('@')
1190        if depvars[0] == 'deaccum': diag='deaccum'
[649]1191        if depvars[0] == 'accum': diag='accum'
[365]1192        for depv in depvars:
1193            if not ncobj.variables.has_key(depv) and not                             \
[756]1194              gen.searchInlist(NONcheckingvars, depv) and                            \
[865]1195              not gen.searchInlist(methods, depv) and not depvars[0] == 'deaccum'    \
1196              and not depvars[0] == 'accum':
[365]1197                print errormsg
1198                print '  ' + main + ": file '" + opts.ncfile +                       \
1199                  "' does not have variable '" + depv + "' !!"
1200                quit(-1)
1201    else:
1202        depvars = diags[idiag].split('|')[1]
1203        if not ncobj.variables.has_key(depvars) and not                              \
[756]1204          gen.searchInlist(NONcheckingvars, depvars) and                             \
1205          not gen.searchInlist(methods, depvars):
1206            print errormsg
[365]1207            print '  ' + main + ": file '" + opts.ncfile +                           \
1208              "' does not have variable '" + depvars + "' !!"
1209            quit(-1)
1210
1211    print "\n    Computing '" + diag + "' from: ", depvars, '...'
1212
1213# acraintot: accumulated total precipitation from WRF RAINC, RAINNC
1214    if diag == 'ACRAINTOT':
1215           
1216        var0 = ncobj.variables[depvars[0]]
1217        var1 = ncobj.variables[depvars[1]]
1218        diagout = var0[:] + var1[:]
1219
1220        dnamesvar = var0.dimensions
1221        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1222
[649]1223        ncvar.insert_variable(ncobj, 'pracc', diagout, dnamesvar, dvnamesvar, newnc)
[365]1224
[649]1225# accum: acumulation of any variable as (Variable, time [as [tunits]
1226#   from/since ....], newvarname)
1227    elif diag == 'accum':
1228
1229        var0 = ncobj.variables[depvars[0]]
1230        var1 = ncobj.variables[depvars[1]]
1231
1232        dnamesvar = var0.dimensions
1233        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1234
1235        diagout, diagoutd, diagoutvd = compute_accum(var0,dnamesvar,dvnamesvar)
1236
1237        CFvarn = ncvar.variables_values(depvars[0])[0]
1238
1239# Removing the flux
1240        if depvars[1] == 'XTIME':
1241            dtimeunits = var1.getncattr('description')
1242            tunits = dtimeunits.split(' ')[0]
1243        else:
1244            dtimeunits = var1.getncattr('units')
1245            tunits = dtimeunits.split(' ')[0]
1246
1247        dtime = (var1[1] - var1[0])*timeunits_seconds(tunits)
1248
1249        ncvar.insert_variable(ncobj, CFvarn + 'acc', diagout*dtime, diagoutd, diagoutvd, newnc)
1250
[365]1251# cllmh with cldfra, pres
1252    elif diag == 'cllmh':
1253           
1254        var0 = ncobj.variables[depvars[0]]
1255        if depvars[1] == 'WRFp':
1256            var1 = WRFp
1257        else:
1258            var01 = ncobj.variables[depvars[1]]
1259            if len(size(var1.shape)) < len(size(var0.shape)):
1260                var1 = np.brodcast_arrays(var01,var0)[0]
1261            else:
1262                var1 = var01
1263
[772]1264        diagout, diagoutd, diagoutvd = Forcompute_cllmh(var0,var1,dnames,dvnames)
1265
[1351]1266        # Removing the nonChecking variable-dimensions from the initial list
1267        varsadd = []
1268        for nonvd in NONchkvardims:
1269            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
1270            varsadd.append(nonvd)
1271
[365]1272        ncvar.insert_variable(ncobj, 'cll', diagout[0,:], diagoutd, diagoutvd, newnc)
1273        ncvar.insert_variable(ncobj, 'clm', diagout[1,:], diagoutd, diagoutvd, newnc)
1274        ncvar.insert_variable(ncobj, 'clh', diagout[2,:], diagoutd, diagoutvd, newnc)
1275
1276# clt with cldfra
1277    elif diag == 'clt':
1278           
1279        var0 = ncobj.variables[depvars]
[772]1280        diagout, diagoutd, diagoutvd = Forcompute_clt(var0,dnames,dvnames)
[1351]1281
1282        # Removing the nonChecking variable-dimensions from the initial list
1283        varsadd = []
1284        for nonvd in NONchkvardims:
1285            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
1286            varsadd.append(nonvd)
1287           
[365]1288        ncvar.insert_variable(ncobj, 'clt', diagout, diagoutd, diagoutvd, newnc)
1289
1290# deaccum: deacumulation of any variable as (Variable, time [as [tunits]
1291#   from/since ....], newvarname)
1292    elif diag == 'deaccum':
1293
1294        var0 = ncobj.variables[depvars[1]]
1295        var1 = ncobj.variables[depvars[2]]
1296
1297        dnamesvar = var0.dimensions
1298        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1299
1300        diagout, diagoutd, diagoutvd = compute_deaccum(var0,dnamesvar,dvnamesvar)
1301
1302# Transforming to a flux
1303        if depvars[2] == 'XTIME':
1304            dtimeunits = var1.getncattr('description')
1305            tunits = dtimeunits.split(' ')[0]
1306        else:
1307            dtimeunits = var1.getncattr('units')
1308            tunits = dtimeunits.split(' ')[0]
1309
1310        dtime = (var1[1] - var1[0])*timeunits_seconds(tunits)
1311        ncvar.insert_variable(ncobj, depvars[3], diagout/dtime, diagoutd, diagoutvd, newnc)
1312
1313# LMDZrh (pres, t, r)
1314    elif diag == 'LMDZrh':
1315           
1316        var0 = ncobj.variables[depvars[0]][:]
1317        var1 = ncobj.variables[depvars[1]][:]
1318        var2 = ncobj.variables[depvars[2]][:]
1319
1320        diagout, diagoutd, diagoutvd = compute_rh(var0,var1,var2,dnames,dvnames)
[1079]1321        ncvar.insert_variable(ncobj, 'hur', diagout, diagoutd, diagoutvd, newnc)
[365]1322
1323# LMDZrhs (psol, t2m, q2m)
1324    elif diag == 'LMDZrhs':
1325           
1326        var0 = ncobj.variables[depvars[0]][:]
1327        var1 = ncobj.variables[depvars[1]][:]
1328        var2 = ncobj.variables[depvars[2]][:]
1329
1330        dnamesvar = ncobj.variables[depvars[0]].dimensions
1331        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1332
1333        diagout, diagoutd, diagoutvd = compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
1334
[1079]1335        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
[365]1336
1337# mslp: mean sea level pressure (pres, psfc, terrain, temp, qv)
1338    elif diag == 'mslp' or diag == 'WRFmslp':
1339           
1340        var1 = ncobj.variables[depvars[1]][:]
1341        var2 = ncobj.variables[depvars[2]][:]
1342        var4 = ncobj.variables[depvars[4]][:]
1343
1344        if diag == 'WRFmslp':
1345            var0 = WRFp
1346            var3 = WRFt
1347            dnamesvar = ncobj.variables['P'].dimensions
1348        else:
1349            var0 = ncobj.variables[depvars[0]][:]
1350            var3 = ncobj.variables[depvars[3]][:]
1351            dnamesvar = ncobj.variables[depvars[0]].dimensions
1352
1353        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1354
1355        diagout, diagoutd, diagoutvd = compute_mslp(var0, var1, var2, var3, var4,    \
1356          dnamesvar, dvnamesvar)
1357
[1581]1358        # Removing the nonChecking variable-dimensions from the initial list
1359        varsadd = []
1360        diagoutvd = list(dvnames)
1361        for nonvd in NONchkvardims:
1362            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1363            varsadd.append(nonvd)
[365]1364        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
1365
[642]1366# OMEGAw (omega, p, t) from NCL formulation (https://www.ncl.ucar.edu/Document/Functions/Contributed/omega_to_w.shtml)
1367    elif diag == 'OMEGAw':
1368           
1369        var0 = ncobj.variables[depvars[0]][:]
1370        var1 = ncobj.variables[depvars[1]][:]
[643]1371        var2 = ncobj.variables[depvars[2]][:]
[642]1372
1373        dnamesvar = ncobj.variables[depvars[0]].dimensions
1374        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1375
1376        diagout, diagoutd, diagoutvd = compute_OMEGAw(var0,var1,var2,dnamesvar,dvnamesvar)
1377
1378        ncvar.insert_variable(ncobj, 'wa', diagout, diagoutd, diagoutvd, newnc)
1379
[365]1380# raintot: instantaneous total precipitation from WRF as (RAINC + RAINC) / dTime
1381    elif diag == 'RAINTOT':
1382
1383        var0 = ncobj.variables[depvars[0]]
1384        var1 = ncobj.variables[depvars[1]]
[445]1385        if depvars[2] != 'WRFtime':
[443]1386            var2 = ncobj.variables[depvars[2]]
[654]1387        else:
1388            var2 = np.arange(var0.shape[0], dtype=int)
[365]1389
1390        var = var0[:] + var1[:]
1391
1392        dnamesvar = var0.dimensions
1393        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1394
1395        diagout, diagoutd, diagoutvd = compute_deaccum(var,dnamesvar,dvnamesvar)
1396
1397# Transforming to a flux
[654]1398        if var2.shape[0] > 1:
[600]1399            if depvars[2] != 'WRFtime':
1400                dtimeunits = var2.getncattr('units')
1401                tunits = dtimeunits.split(' ')[0]
1402   
1403                dtime = (var2[1] - var2[0])*timeunits_seconds(tunits)
1404            else:
1405                var2 = ncobj.variables['Times']
1406                time1 = var2[0,:]
1407                time2 = var2[1,:]
1408                tmf1 = ''
1409                tmf2 = ''
1410                for ic in range(len(time1)):
1411                    tmf1 = tmf1 + time1[ic]
1412                    tmf2 = tmf2 + time2[ic]
[654]1413                dtdate1 = dtime.datetime.strptime(tmf1,"%Y-%m-%d_%H:%M:%S")
1414                dtdate2 = dtime.datetime.strptime(tmf2,"%Y-%m-%d_%H:%M:%S")
[600]1415                diffdate12 = dtdate2 - dtdate1
1416                dtime = diffdate12.total_seconds()
1417                print 'dtime:',dtime
[442]1418        else:
[600]1419            print warnmsg
1420            print '  ' + fname + ": only 1 time-step for '" + diag + "' !!"
1421            print '    leaving a zero value!'
1422            diagout = var0*0.
1423            dtime=1.
[442]1424
[365]1425        ncvar.insert_variable(ncobj, 'pr', diagout/dtime, diagoutd, diagoutvd, newnc)
1426
[612]1427# rhs (psfc, t, q) from TimeSeries files
1428    elif diag == 'TSrhs':
1429           
1430        p0=100000.
1431        var0 = ncobj.variables[depvars[0]][:]
1432        var1 = (ncobj.variables[depvars[1]][:])*(var0/p0)**(2./7.)
1433        var2 = ncobj.variables[depvars[2]][:]
1434
1435        dnamesvar = ncobj.variables[depvars[0]].dimensions
1436        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1437
1438        diagout, diagoutd, diagoutvd = compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
1439
[1079]1440        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
[612]1441
1442# td (psfc, t, q) from TimeSeries files
[613]1443    elif diag == 'TStd' or diag == 'td':
[612]1444           
1445        var0 = ncobj.variables[depvars[0]][:]
1446        var1 = ncobj.variables[depvars[1]][:] - 273.15
1447        var2 = ncobj.variables[depvars[2]][:]
1448
1449        dnamesvar = ncobj.variables[depvars[0]].dimensions
1450        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1451
1452        diagout, diagoutd, diagoutvd = compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
1453
1454        ncvar.insert_variable(ncobj, 'tds', diagout, diagoutd, diagoutvd, newnc)
1455
1456# td (psfc, t, q) from TimeSeries files
[616]1457    elif diag == 'TStdC' or diag == 'tdC':
[612]1458           
1459        var0 = ncobj.variables[depvars[0]][:]
1460# Temperature is already in degrees Celsius
1461        var1 = ncobj.variables[depvars[1]][:]
1462        var2 = ncobj.variables[depvars[2]][:]
1463
1464        dnamesvar = ncobj.variables[depvars[0]].dimensions
1465        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1466
1467        diagout, diagoutd, diagoutvd = compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
1468
1469        ncvar.insert_variable(ncobj, 'tds', diagout, diagoutd, diagoutvd, newnc)
1470
1471# wds (u, v)
1472    elif diag == 'TSwds' or diag == 'wds' :
1473 
1474        var0 = ncobj.variables[depvars[0]][:]
1475        var1 = ncobj.variables[depvars[1]][:]
1476
1477        dnamesvar = ncobj.variables[depvars[0]].dimensions
1478        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1479
1480        diagout, diagoutd, diagoutvd = compute_wds(var0,var1,dnamesvar,dvnamesvar)
1481
1482        ncvar.insert_variable(ncobj, 'wds', diagout, diagoutd, diagoutvd, newnc)
1483
1484# wss (u, v)
[613]1485    elif diag == 'TSwss' or diag == 'wss':
[612]1486           
1487        var0 = ncobj.variables[depvars[0]][:]
1488        var1 = ncobj.variables[depvars[1]][:]
1489
1490        dnamesvar = ncobj.variables[depvars[0]].dimensions
1491        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1492
1493        diagout, diagoutd, diagoutvd = compute_wss(var0,var1,dnamesvar,dvnamesvar)
1494
1495        ncvar.insert_variable(ncobj, 'wss', diagout, diagoutd, diagoutvd, newnc)
1496
[365]1497# turbulence (var)
1498    elif diag == 'turbulence':
1499
1500        var0 = ncobj.variables[depvars][:]
1501
1502        dnamesvar = list(ncobj.variables[depvars].dimensions)
1503        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1504
1505        diagout, diagoutd, diagoutvd = compute_turbulence(var0,dnamesvar,dvnamesvar)
[959]1506        valsvar = gen.variables_values(depvars)
[365]1507
[959]1508        newvarn = depvars + 'turb'
1509        print main + '; Lluis newvarn:', newvarn
1510        ncvar.insert_variable(ncobj, newvarn, diagout, diagoutd, 
[365]1511          diagoutvd, newnc)
[959]1512        print main + '; Lluis variables:', newnc.variables.keys()
1513        varobj = newnc.variables[newvarn]
[365]1514        attrv = varobj.long_name
1515        attr = varobj.delncattr('long_name')
1516        newattr = ncvar.set_attribute(varobj, 'long_name', attrv +                   \
1517          " Taylor decomposition turbulence term")
1518
[390]1519# WRFbils fom WRF as HFX + LH
1520    elif diag == 'WRFbils':
1521           
1522        var0 = ncobj.variables[depvars[0]][:]
1523        var1 = ncobj.variables[depvars[1]][:]
1524
1525        diagout = var0 + var1
[867]1526        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1527        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
[390]1528
[867]1529        ncvar.insert_variable(ncobj, 'bils', diagout, dnamesvar, dvnamesvar, newnc)
[390]1530
[1581]1531# WRFclivi WRF water vapour path WRFdens, QICE, QGRAUPEL, QHAIL
1532    elif diag == 'WRFclivi':
1533           
1534        var0 = WRFdens
1535        qtot = ncobj.variables[depvars[1]]
1536        qtotv = qtot[:]
1537        Nspecies = len(depvars) - 2
1538        for iv in range(Nspecies):
[1585]1539            if ncobj.variables.has_key(depvars[iv+2]):
1540                var1 = ncobj.variables[depvars[iv+2]][:]
1541                qtotv = qtotv + var1
[1581]1542
1543        dnamesvar = list(qtot.dimensions)
1544        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1545
1546        diagout, diagoutd, diagoutvd = compute_clivi(var0, qtotv, dnamesvar,dvnamesvar)
1547
1548        # Removing the nonChecking variable-dimensions from the initial list
1549        varsadd = []
1550        diagoutvd = list(dvnames)
1551        for nonvd in NONchkvardims:
1552            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1553            varsadd.append(nonvd)
1554        ncvar.insert_variable(ncobj, 'clivi', diagout, diagoutd, diagoutvd, newnc)
1555
1556# WRFclwvl WRF water cloud-condensed path WRFdens, QCLOUD, QICE, QGRAUPEL, QHAIL
1557    elif diag == 'WRFclwvl':
1558           
1559        var0 = WRFdens
1560        qtot = ncobj.variables[depvars[1]]
1561        qtotv = ncobj.variables[depvars[1]]
1562        Nspecies = len(depvars) - 2
1563        for iv in range(Nspecies):
[1585]1564            if ncobj.variables.has_key(depvars[iv+2]):
1565                var1 = ncobj.variables[depvars[iv+2]]
1566                qtotv = qtotv + var1[:]
[1581]1567
1568        dnamesvar = list(qtot.dimensions)
1569        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1570
1571        diagout, diagoutd, diagoutvd = compute_clwvl(var0, qtotv, dnamesvar,dvnamesvar)
1572
1573        # Removing the nonChecking variable-dimensions from the initial list
1574        varsadd = []
1575        diagoutvd = list(dvnames)
1576        for nonvd in NONchkvardims:
1577            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1578            varsadd.append(nonvd)
1579        ncvar.insert_variable(ncobj, 'clwvl', diagout, diagoutd, diagoutvd, newnc)
1580
[654]1581# WRFgeop geopotential from WRF as PH + PHB
1582    elif diag == 'WRFgeop':
[1382]1583        var0 = ncobj.variables[depvars[0]][:]
1584        var1 = ncobj.variables[depvars[1]][:]
[654]1585
[1382]1586        # de-staggering geopotential
1587        diagout0 = var0 + var1
1588        dt = diagout0.shape[0]
1589        dz = diagout0.shape[1]
1590        dy = diagout0.shape[2]
1591        dx = diagout0.shape[3]
1592
1593        diagout = np.zeros((dt,dz-1,dy,dx), dtype=np.float)
1594        diagout = 0.5*(diagout0[:,1:dz,:,:]+diagout0[:,0:dz-1,:,:])
1595
1596        # Removing the nonChecking variable-dimensions from the initial list
1597        varsadd = []
[1389]1598        diagoutvd = list(dvnames)
[1382]1599        for nonvd in NONchkvardims:
[1389]1600            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
[1382]1601            varsadd.append(nonvd)
1602
[1389]1603        ncvar.insert_variable(ncobj, 'zg', diagout, dnames, diagoutvd, newnc)
[654]1604
[390]1605# WRFp pressure from WRF as P + PB
[365]1606    elif diag == 'WRFp':
1607           
1608        diagout = WRFp
1609
1610        ncvar.insert_variable(ncobj, 'pres', diagout, dnames, dvnames, newnc)
1611
1612# WRFpos
1613    elif diag == 'WRFpos':
1614           
1615        dnamesvar = ncobj.variables['MAPFAC_M'].dimensions
1616        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1617
1618        ncvar.insert_variable(ncobj, 'WRFpos', WRFpos, dnamesvar, dvnamesvar, newnc)
1619
1620# WRFprw WRF water vapour path WRFdens, QVAPOR
1621    elif diag == 'WRFprw':
1622           
1623        var0 = WRFdens
1624        var1 = ncobj.variables[depvars[1]]
1625
1626        dnamesvar = list(var1.dimensions)
1627        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1628
1629        diagout, diagoutd, diagoutvd = compute_prw(var0, var1, dnamesvar,dvnamesvar)
1630
[1586]1631        # Removing the nonChecking variable-dimensions from the initial list
1632        varsadd = []
1633        diagoutvd = list(dvnames)
1634        for nonvd in NONchkvardims:
1635            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1636            varsadd.append(nonvd)
[365]1637        ncvar.insert_variable(ncobj, 'prw', diagout, diagoutd, diagoutvd, newnc)
1638
1639# WRFrh (P, T, QVAPOR)
1640    elif diag == 'WRFrh':
1641           
1642        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
1643        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1644
[878]1645        ncvar.insert_variable(ncobj, 'hur', WRFrh, dnames, dvnames, newnc)
[365]1646
1647# WRFrhs (PSFC, T2, Q2)
1648    elif diag == 'WRFrhs':
1649           
1650        var0 = ncobj.variables[depvars[0]][:]
1651        var1 = ncobj.variables[depvars[1]][:]
1652        var2 = ncobj.variables[depvars[2]][:]
1653
1654        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
1655        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1656
1657        diagout, diagoutd, diagoutvd = compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
[878]1658        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
[365]1659
1660# rvors (u10, v10, WRFpos)
1661    elif diag == 'WRFrvors':
1662           
1663        var0 = ncobj.variables[depvars[0]]
1664        var1 = ncobj.variables[depvars[1]]
1665
1666        diagout = rotational_z(var0, var1, distx)
1667
1668        dnamesvar = ncobj.variables[depvars[0]].dimensions
1669        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1670
1671        ncvar.insert_variable(ncobj, 'rvors', diagout, dnamesvar, dvnamesvar, newnc)
1672
[884]1673# WRFt (T, P, PB)
1674    elif diag == 'WRFt':
1675        var0 = ncobj.variables[depvars[0]][:]
1676        var1 = ncobj.variables[depvars[1]][:]
1677        var2 = ncobj.variables[depvars[2]][:]
[654]1678
[884]1679        p0=100000.
1680        p=var1 + var2
1681
1682        WRFt = (var0 + 300.)*(p/p0)**(2./7.)
1683
1684        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
1685        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1686
[1382]1687        # Removing the nonChecking variable-dimensions from the initial list
1688        varsadd = []
[1389]1689        diagoutvd = list(dvnames)
[1382]1690        for nonvd in NONchkvardims:
[1389]1691            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
[1382]1692            varsadd.append(nonvd)
1693
[1389]1694        ncvar.insert_variable(ncobj, 'ta', WRFt, dnames, diagoutvd, newnc)
[884]1695
[914]1696# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
1697    elif diag == 'WRFua':
1698        var0 = ncobj.variables[depvars[0]][:]
1699        var1 = ncobj.variables[depvars[1]][:]
1700        var2 = ncobj.variables[depvars[2]][:]
1701        var3 = ncobj.variables[depvars[3]][:]
1702
1703        # un-staggering variables
1704        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1705        ua = np.zeros(tuple(unstgdims), dtype=np.float)
1706        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1707        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1708        unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1709        unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1710
1711        for iz in range(var0.shape[1]):
1712            ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
1713
[959]1714        dnamesvar = ['Time','bottom_top','south_north','west_east']
[914]1715        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1716
[1404]1717        # Removing the nonChecking variable-dimensions from the initial list
1718        varsadd = []
1719        diagoutvd = list(dvnames)
1720        for nonvd in NONchkvardims:
1721            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1722            varsadd.append(nonvd)
[914]1723
[1404]1724        ncvar.insert_variable(ncobj, 'ua', ua, dnames, diagoutvd, newnc)
1725
[914]1726# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
1727    elif diag == 'WRFva':
1728        var0 = ncobj.variables[depvars[0]][:]
1729        var1 = ncobj.variables[depvars[1]][:]
1730        var2 = ncobj.variables[depvars[2]][:]
1731        var3 = ncobj.variables[depvars[3]][:]
1732
1733        # un-staggering variables
1734        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1735        va = np.zeros(tuple(unstgdims), dtype=np.float)
1736        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1737        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1738        unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1739        unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1740        for iz in range(var0.shape[1]):
1741            va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
1742
[959]1743        dnamesvar = ['Time','bottom_top','south_north','west_east']
[914]1744        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1745
[1404]1746        # Removing the nonChecking variable-dimensions from the initial list
1747        varsadd = []
1748        diagoutvd = list(dvnames)
1749        for nonvd in NONchkvardims:
1750            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1751            varsadd.append(nonvd)
1752        ncvar.insert_variable(ncobj, 'va', va, dnames, diagoutvd, newnc)
[914]1753
1754# WRFtime
[654]1755    elif diag == 'WRFtime':
1756           
1757        diagout = WRFtime
1758
1759        dnamesvar = ['Time']
1760        dvnamesvar = ['Times']
1761
1762        ncvar.insert_variable(ncobj, 'time', diagout, dnamesvar, dvnamesvar, newnc)
1763
[959]1764# ws (U, V)
1765    elif diag == 'ws':
1766           
1767        var0 = ncobj.variables[depvars[0]][:]
1768        var1 = ncobj.variables[depvars[1]][:]
1769        # un-staggering variables
1770        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1771        va = np.zeros(tuple(unstgdims), dtype=np.float)
1772        unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1773        unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1774        unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1775        unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1776
1777        dnamesvar = ['Time','bottom_top','south_north','west_east']
1778        diagout = np.sqrt(unstgvar0*unstgvar0 + unstgvar1*unstgvar1)
1779
1780#        dnamesvar = ncobj.variables[depvars[0]].dimensions
1781        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1782
[1408]1783        # Removing the nonChecking variable-dimensions from the initial list
1784        varsadd = []
1785        diagoutvd = list(dvnamesvar)
1786        for nonvd in NONchkvardims:
1787            if gen.searchInlist(dvnamesvar,nonvd): diagoutvd.remove(nonvd)
1788            varsadd.append(nonvd)
1789        ncvar.insert_variable(ncobj, 'ws', diagout, dnamesvar, diagoutvd, newnc)
[959]1790
[365]1791# wss (u10, v10)
1792    elif diag == 'wss':
1793           
1794        var0 = ncobj.variables[depvars[0]][:]
1795        var1 = ncobj.variables[depvars[1]][:]
1796
1797        diagout = np.sqrt(var0*var0 + var1*var1)
1798
1799        dnamesvar = ncobj.variables[depvars[0]].dimensions
1800        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
1801
1802        ncvar.insert_variable(ncobj, 'wss', diagout, dnamesvar, dvnamesvar, newnc)
1803
[654]1804# WRFheight height from WRF geopotential as WRFGeop/g
1805    elif diag == 'WRFheight':
1806           
1807        diagout = WRFgeop/grav
1808
[1412]1809        # Removing the nonChecking variable-dimensions from the initial list
1810        varsadd = []
1811        diagoutvd = list(dvnames)
1812        for nonvd in NONchkvardims:
1813            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1814            varsadd.append(nonvd)
[654]1815
[1412]1816        ncvar.insert_variable(ncobj, 'zhgt', diagout, dnames, diagoutvd, newnc)
1817
[1413]1818# WRFheightrel relative-height from WRF geopotential as WRFgeop(PH + PHB)/g-HGT 'WRFheightrel|PH@PHB@HGT
1819    elif diag == 'WRFheightrel':
1820        var0 = ncobj.variables[depvars[0]][:]
1821        var1 = ncobj.variables[depvars[1]][:]
1822        var2 = ncobj.variables[depvars[2]][:]
1823
1824        dimz = var0.shape[1]
1825        diagout = np.zeros(tuple(var0.shape), dtype=np.float)
1826        for iz in range(dimz):
[1419]1827            diagout[:,iz,:,:] = (var0[:,iz,:,:]+ var1[:,iz,:,:])/grav - var2
[1413]1828
1829        # Removing the nonChecking variable-dimensions from the initial list
1830        varsadd = []
1831        diagoutvd = list(dvnames)
1832        for nonvd in NONchkvardims:
1833            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
1834            varsadd.append(nonvd)
1835
1836        ncvar.insert_variable(ncobj, 'zhgtrel', diagout, dnames, diagoutvd, newnc)
1837
[365]1838    else:
1839        print errormsg
1840        print '  ' + main + ": diagnostic '" + diag + "' not ready!!!"
1841        print '    available diagnostics: ', availdiags
1842        quit(-1)
1843
1844    newnc.sync()
[1351]1845    # Adding that additional variables required to compute some diagnostics which
1846    #   where not in the original file
1847    for vadd in varsadd:
1848        if not gen.searchInlist(newnc.variables.keys(),vadd):
1849            attrs = dictcompvars[vadd]
1850            vvn = attrs['name']
1851            if not gen.searchInlist(newnc.variables.keys(), vvn):
1852                iidvn = dvnames.index(vadd)
1853                dnn = dnames[iidvn]
1854                if vadd == 'WRFtime':
1855                    dvarvals = WRFtime[:]
1856                newvar = newnc.createVariable(vvn, 'f8', (dnn))
1857                newvar[:] = dvarvals
1858                for attn in attrs.keys():
1859                    if attn != 'name':
1860                        attv = attrs[attn]
1861                        ncvar.set_attribute(newvar, attn, attv)
[365]1862
1863#   end of diagnostics
1864
1865# Global attributes
1866##
1867atvar = ncvar.set_attribute(newnc, 'program', 'diagnostics.py')
1868atvar = ncvar.set_attribute(newnc, 'version', '1.0')
1869atvar = ncvar.set_attribute(newnc, 'author', 'Fita Borrell, Lluis')
1870atvar = ncvar.set_attribute(newnc, 'institution', 'Laboratoire Meteorologie ' +      \
1871  'Dynamique')
1872atvar = ncvar.set_attribute(newnc, 'university', 'Universite Pierre et Marie ' +     \
1873  'Curie -- Jussieu')
1874atvar = ncvar.set_attribute(newnc, 'centre', 'Centre national de la recherche ' +    \
1875  'scientifique')
1876atvar = ncvar.set_attribute(newnc, 'city', 'Paris')
1877atvar = ncvar.set_attribute(newnc, 'original_file', opts.ncfile)
1878
1879gorigattrs = ncobj.ncattrs()
1880
1881for attr in gorigattrs:
1882    attrv = ncobj.getncattr(attr)
1883    atvar = ncvar.set_attribute(newnc, attr, attrv)
1884
1885ncobj.close()
1886newnc.close()
1887
1888print '\n' + main + ': successfull writting of diagnostics file "' + ofile + '" !!!'
Note: See TracBrowser for help on using the repository browser.