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

Last change on this file since 653 was 649, checked in by lfita, 10 years ago

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