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

Last change on this file since 664 was 654, checked in by lfita, 9 years ago

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