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

Last change on this file since 762 was 756, checked in by lfita, 9 years ago

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