source: lmdz_wrf/trunk/tools/diag_tools.py @ 1700

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

Adding filling of 4D pressure variable for computation of 'td'
File size: 54.9 KB
Line 
1# Tools for the compute of diagnostics
2# L. Fita, CIMA. CONICET-UBA, CNRS UMI-IFAECI, Buenos Aires, Argentina
3
4# Available general pupose diagnostics (model independent) providing (varv1, varv2, ..., dimns, dimvns)
5# compute_accum: Function to compute the accumulation of a variable
6# compute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following
7#   newmicro.F90 from LMDZ compute_clt(cldfra, pres, dimns, dimvns)
8# compute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ
9# compute_clivi: Function to compute cloud-ice water path (clivi)
10# compute_clwvl: Function to compute condensed water path (clwvl)
11# compute_deaccum: Function to compute the deaccumulation of a variable
12# compute_mslp: Function to compute mslp: mean sea level pressure following p_interp.F90 from WRF
13# compute_OMEGAw: Function to transform OMEGA [Pas-1] to velocities [ms-1]
14# compute_prw: Function to compute water vapour path (prw)
15# compute_rh: Function to compute relative humidity following 'Tetens' equation (T,P) ...'
16# compute_td: Function to compute the dew point temperature
17# compute_turbulence: Function to compute the rubulence term of the Taylor's decomposition ...'
18# C_diagnostic: Class to compute generic variables
19# compute_wds: Function to compute the wind direction
20# compute_wss: Function to compute the wind speed
21# compute_WRFta: Function to compute WRF air temperature
22# compute_WRFtd: Function to compute WRF dew-point air temperature
23# compute_WRFua: Function to compute geographical rotated WRF x-wind
24# compute_WRFva: Function to compute geographical rotated WRF y-wind
25# compute_WRFuava: Function to compute geographical rotated WRF 3D winds
26# compute_WRFuas: Function to compute geographical rotated WRF 2-meter x-wind
27# compute_WRFvas: Function to compute geographical rotated WRF 2-meter y-wind
28# compute_WRFuasvas: Fucntion to compute geographical rotated WRF 2-meter winds
29# derivate_centered: Function to compute the centered derivate of a given field
30# def Forcompute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ via Fortran subroutine
31# Forcompute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ via a Fortran module
32# W_diagnostic: Class to compute WRF diagnostics variables
33
34# Others just providing variable values
35# var_cllmh: Fcuntion to compute cllmh on a 1D column
36# var_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ using 1D vertical column values
37# var_mslp: Fcuntion to compute mean sea-level pressure
38# var_td: Function to compute dew-point air temperature from temperature and pressure values
39# var_virtualTemp: This function returns virtual temperature in K,
40# var_WRFtime: Function to copmute CFtimes from WRFtime variable
41# var_wd: Function to compute the wind direction
42# var_wd: Function to compute the wind speed
43# rotational_z: z-component of the rotatinoal of horizontal vectorial field
44# turbulence_var: Function to compute the Taylor's decomposition turbulence term from a a given variable
45
46import numpy as np
47from netCDF4 import Dataset as NetCDFFile
48import os
49import re
50import nc_var_tools as ncvar
51import generic_tools as gen
52import datetime as dtime
53import module_ForDiag as fdin
54import module_ForDef as fdef
55
56main = 'diag_tools.py'
57errormsg = 'ERROR -- error -- ERROR -- error'
58warnmsg = 'WARNING -- warning -- WARNING -- warning'
59
60# Constants
61grav = fdef.module_definitions.grav
62
63# Available WRFiag
64Wavailablediags = ['p', 'ta', 'td', 'ua', 'va', 'uas', 'vas', 'wd', 'ws', 'zg']
65
66# Available General diagnostics
67Cavailablediags = ['td', 'wd', 'ws']
68
69# Gneral information
70##
71def reduce_spaces(string):
72    """ Function to give words of a line of text removing any extra space
73    """
74    values = string.replace('\n','').split(' ')
75    vals = []
76    for val in values:
77         if len(val) > 0:
78             vals.append(val)
79
80    return vals
81
82def variable_combo(varn,combofile):
83    """ Function to provide variables combination from a given variable name
84      varn= name of the variable
85      combofile= ASCII file with the combination of variables
86        [varn] [combo]
87          [combo]: '@' separated list of variables to use to generate [varn]
88            [WRFdt] to get WRF time-step (from general attributes)
89    >>> variable_combo('WRFprls','/home/lluis/PY/diagnostics.inf')
90    deaccum@RAINNC@XTIME@prnc
91    """
92    fname = 'variable_combo'
93
94    if varn == 'h':
95        print fname + '_____________________________________________________________'
96        print variable_combo.__doc__
97        quit()
98
99    if not os.path.isfile(combofile):
100        print errormsg
101        print '  ' + fname + ": file with combinations '" + combofile +              \
102          "' does not exist!!"
103        quit(-1)
104
105    objf = open(combofile, 'r')
106
107    found = False
108    for line in objf:
109        linevals = reduce_spaces(line)
110        varnf = linevals[0]
111        combo = linevals[1].replace('\n','')
112        if varn == varnf: 
113            found = True
114            break
115
116    if not found:
117        print errormsg
118        print '  ' + fname + ": variable '" + varn + "' not found in '" + combofile +\
119          "' !!"
120        combo='ERROR'
121
122    objf.close()
123
124    return combo
125
126# Mathematical operators
127##
128def compute_accum(varv, dimns, dimvns):
129    """ Function to compute the accumulation of a variable
130    compute_accum(varv, dimnames, dimvns)
131      [varv]= values to accum (assuming [t,])
132      [dimns]= list of the name of the dimensions of the [varv]
133      [dimvns]= list of the name of the variables with the values of the
134        dimensions of [varv]
135    """
136    fname = 'compute_accum'
137
138    deacdims = dimns[:]
139    deacvdims = dimvns[:]
140
141    slicei = []
142    slicee = []
143
144    Ndims = len(varv.shape)
145    for iid in range(0,Ndims):
146        slicei.append(slice(0,varv.shape[iid]))
147        slicee.append(slice(0,varv.shape[iid]))
148
149    slicee[0] = np.arange(varv.shape[0])
150    slicei[0] = np.arange(varv.shape[0])
151    slicei[0][1:varv.shape[0]] = np.arange(varv.shape[0]-1)
152
153    vari = varv[tuple(slicei)]
154    vare = varv[tuple(slicee)]
155
156    ac = vari*0.
157    for it in range(1,varv.shape[0]):
158        ac[it,] = ac[it-1,] + vare[it,]
159
160    return ac, deacdims, deacvdims
161
162def compute_deaccum(varv, dimns, dimvns):
163    """ Function to compute the deaccumulation of a variable
164    compute_deaccum(varv, dimnames, dimvns)
165      [varv]= values to deaccum (assuming [t,])
166      [dimns]= list of the name of the dimensions of the [varv]
167      [dimvns]= list of the name of the variables with the values of the
168        dimensions of [varv]
169    """
170    fname = 'compute_deaccum'
171
172    deacdims = dimns[:]
173    deacvdims = dimvns[:]
174
175    slicei = []
176    slicee = []
177
178    Ndims = len(varv.shape)
179    for iid in range(0,Ndims):
180        slicei.append(slice(0,varv.shape[iid]))
181        slicee.append(slice(0,varv.shape[iid]))
182
183    slicee[0] = np.arange(varv.shape[0])
184    slicei[0] = np.arange(varv.shape[0])
185    slicei[0][1:varv.shape[0]] = np.arange(varv.shape[0]-1)
186
187    vari = varv[tuple(slicei)]
188    vare = varv[tuple(slicee)]
189
190    deac = vare - vari
191
192    return deac, deacdims, deacvdims
193
194def derivate_centered(var,dim,dimv):
195    """ Function to compute the centered derivate of a given field
196      centered derivate(n) = (var(n-1) + var(n+1))/(2*dn).
197    [var]= variable
198    [dim]= which dimension to compute the derivate
199    [dimv]= dimension values (can be of different dimension of [var])
200    >>> derivate_centered(np.arange(16).reshape(4,4)*1.,1,1.)
201    [[  0.   1.   2.   0.]
202     [  0.   5.   6.   0.]
203     [  0.   9.  10.   0.]
204     [  0.  13.  14.   0.]]
205    """
206
207    fname = 'derivate_centered'
208   
209    vark = var.dtype
210
211    if hasattr(dimv, "__len__"):
212# Assuming that the last dimensions of var [..., N, M] are the same of dimv [N, M]
213        if len(var.shape) != len(dimv.shape):
214            dimvals = np.zeros((var.shape), dtype=vark)
215            if len(var.shape) - len(dimv.shape) == 1:
216                for iz in range(var.shape[0]):
217                    dimvals[iz,] = dimv
218            elif len(var.shape) - len(dimv.shape) == 2:
219                for it in range(var.shape[0]):
220                    for iz in range(var.shape[1]):
221                        dimvals[it,iz,] = dimv
222            else:
223                print errormsg
224                print '  ' + fname + ': dimension difference between variable',      \
225                  var.shape,'and variable with dimension values',dimv.shape,         \
226                  ' not ready !!!'
227                quit(-1)
228        else:
229            dimvals = dimv
230    else:
231# dimension values are identical everywhere!
232# from: http://stackoverflow.com/questions/16807011/python-how-to-identify-if-a-variable-is-an-array-or-a-scalar   
233        dimvals = np.ones((var.shape), dtype=vark)*dimv
234
235    derivate = np.zeros((var.shape), dtype=vark)
236    if dim > len(var.shape) - 1:
237        print errormsg
238        print '  ' + fname + ': dimension',dim,' too big for given variable of ' +   \
239          'shape:', var.shape,'!!!'
240        quit(-1)
241
242    slicebef = []
243    sliceaft = []
244    sliceder = []
245
246    for id in range(len(var.shape)):
247        if id == dim:
248            slicebef.append(slice(0,var.shape[id]-2))
249            sliceaft.append(slice(2,var.shape[id]))
250            sliceder.append(slice(1,var.shape[id]-1))
251        else:
252            slicebef.append(slice(0,var.shape[id]))
253            sliceaft.append(slice(0,var.shape[id]))
254            sliceder.append(slice(0,var.shape[id]))
255
256    if hasattr(dimv, "__len__"):
257        derivate[tuple(sliceder)] = (var[tuple(slicebef)] + var[tuple(sliceaft)])/   \
258          ((dimvals[tuple(sliceaft)] - dimvals[tuple(slicebef)]))
259        print (dimvals[tuple(sliceaft)] - dimvals[tuple(slicebef)])
260    else:
261        derivate[tuple(sliceder)] = (var[tuple(slicebef)] + var[tuple(sliceaft)])/   \
262          (2.*dimv)
263
264#    print 'before________'
265#    print var[tuple(slicebef)]
266
267#    print 'after________'
268#    print var[tuple(sliceaft)]
269
270    return derivate
271
272def rotational_z(Vx,Vy,pos):
273    """ z-component of the rotatinoal of horizontal vectorial field
274    \/ x (Vx,Vy,Vz) = \/xVy - \/yVx
275    [Vx]= Variable component x
276    [Vy]=  Variable component y
277    [pos]= poisition of the grid points
278    >>> rotational_z(np.arange(16).reshape(4,4)*1., np.arange(16).reshape(4,4)*1., 1.)
279    [[  0.   1.   2.   0.]
280     [ -4.   0.   0.  -7.]
281     [ -8.   0.   0. -11.]
282     [  0.  13.  14.   0.]]
283    """
284
285    fname =  'rotational_z'
286
287    ndims = len(Vx.shape)
288    rot1 = derivate_centered(Vy,ndims-1,pos)
289    rot2 = derivate_centered(Vx,ndims-2,pos)
290
291    rot = rot1 - rot2
292
293    return rot
294
295# Diagnostics
296##
297
298def var_clt(cfra):
299    """ Function to compute the total cloud fraction following 'newmicro.F90' from
300      LMDZ using 1D vertical column values
301      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
302    """
303    ZEPSEC=1.0E-12
304
305    fname = 'var_clt'
306
307    zclear = 1.
308    zcloud = 0.
309
310    dz = cfra.shape[0]
311    for iz in range(dz):
312        zclear =zclear*(1.-np.max([cfra[iz],zcloud]))/(1.-np.min([zcloud,1.-ZEPSEC]))
313        clt = 1. - zclear
314        zcloud = cfra[iz]
315
316    return clt
317
318def compute_clt(cldfra, dimns, dimvns):
319    """ Function to compute the total cloud fraction following 'newmicro.F90' from
320      LMDZ
321    compute_clt(cldfra, dimnames)
322      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
323      [dimns]= list of the name of the dimensions of [cldfra]
324      [dimvns]= list of the name of the variables with the values of the
325        dimensions of [cldfra]
326    """
327    fname = 'compute_clt'
328
329    cltdims = dimns[:]
330    cltvdims = dimvns[:]
331
332    if len(cldfra.shape) == 4:
333        clt = np.zeros((cldfra.shape[0],cldfra.shape[2],cldfra.shape[3]),            \
334          dtype=np.float)
335        dx = cldfra.shape[3]
336        dy = cldfra.shape[2]
337        dz = cldfra.shape[1]
338        dt = cldfra.shape[0]
339        cltdims.pop(1)
340        cltvdims.pop(1)
341
342        for it in range(dt):
343            for ix in range(dx):
344                for iy in range(dy):
345                    zclear = 1.
346                    zcloud = 0.
347                    gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
348                    clt[it,iy,ix] = var_clt(cldfra[it,:,iy,ix])
349
350    else:
351        clt = np.zeros((cldfra.shape[1],cldfra.shape[2]), dtype=np.float)
352        dx = cldfra.shape[2]
353        dy = cldfra.shape[1]
354        dy = cldfra.shape[0]
355        cltdims.pop(0)
356        cltvdims.pop(0)
357        for ix in range(dx):
358            for iy in range(dy):
359                zclear = 1.
360                zcloud = 0.
361                gen.percendone(ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
362                clt[iy,ix] = var_clt(cldfra[:,iy,ix])
363
364    return clt, cltdims, cltvdims
365
366def Forcompute_clt(cldfra, dimns, dimvns):
367    """ Function to compute the total cloud fraction following 'newmicro.F90' from
368      LMDZ via a Fortran module
369    compute_clt(cldfra, dimnames)
370      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
371      [dimns]= list of the name of the dimensions of [cldfra]
372      [dimvns]= list of the name of the variables with the values of the
373        dimensions of [cldfra]
374    """
375    fname = 'Forcompute_clt'
376
377    cltdims = dimns[:]
378    cltvdims = dimvns[:]
379
380
381    if len(cldfra.shape) == 4:
382        clt = np.zeros((cldfra.shape[0],cldfra.shape[2],cldfra.shape[3]),            \
383          dtype=np.float)
384        dx = cldfra.shape[3]
385        dy = cldfra.shape[2]
386        dz = cldfra.shape[1]
387        dt = cldfra.shape[0]
388        cltdims.pop(1)
389        cltvdims.pop(1)
390
391        clt = fdin.module_fordiagnostics.compute_clt4d2(cldfra[:])
392
393    else:
394        clt = np.zeros((cldfra.shape[1],cldfra.shape[2]), dtype=np.float)
395        dx = cldfra.shape[2]
396        dy = cldfra.shape[1]
397        dy = cldfra.shape[0]
398        cltdims.pop(0)
399        cltvdims.pop(0)
400
401        clt = fdin.module_fordiagnostics.compute_clt3d1(cldfra[:])
402
403    return clt, cltdims, cltvdims
404
405def var_cllmh(cfra, p):
406    """ Fcuntion to compute cllmh on a 1D column
407    """
408
409    fname = 'var_cllmh'
410
411    ZEPSEC =1.0E-12
412    prmhc = 440.*100.
413    prmlc = 680.*100.
414
415    zclearl = 1.
416    zcloudl = 0.
417    zclearm = 1.
418    zcloudm = 0.
419    zclearh = 1.
420    zcloudh = 0.
421
422    dvz = cfra.shape[0]
423
424    cllmh = np.ones((3), dtype=np.float)
425
426    for iz in range(dvz):
427        if p[iz] < prmhc:
428            cllmh[2] = cllmh[2]*(1.-np.max([cfra[iz], zcloudh]))/(1.-                \
429              np.min([zcloudh,1.-ZEPSEC]))
430            zcloudh = cfra[iz]
431        elif p[iz] >= prmhc and p[iz] < prmlc:
432            cllmh[1] = cllmh[1]*(1.-np.max([cfra[iz], zcloudm]))/(1.-                \
433              np.min([zcloudm,1.-ZEPSEC]))
434            zcloudm = cfra[iz]
435        elif p[iz] >= prmlc:
436            cllmh[0] = cllmh[0]*(1.-np.max([cfra[iz], zcloudl]))/(1.-                \
437              np.min([zcloudl,1.-ZEPSEC]))
438            zcloudl = cfra[iz]
439
440    cllmh = 1.- cllmh
441
442    return cllmh
443
444def Forcompute_cllmh(cldfra, pres, dimns, dimvns):
445    """ Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ via Fortran subroutine
446    compute_clt(cldfra, pres, dimns, dimvns)
447      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
448      [pres] = pressure field
449      [dimns]= list of the name of the dimensions of [cldfra]
450      [dimvns]= list of the name of the variables with the values of the
451        dimensions of [cldfra]
452    """
453    fname = 'Forcompute_cllmh'
454
455    cllmhdims = dimns[:]
456    cllmhvdims = dimvns[:]
457
458    if len(cldfra.shape) == 4:
459        dx = cldfra.shape[3]
460        dy = cldfra.shape[2]
461        dz = cldfra.shape[1]
462        dt = cldfra.shape[0]
463        cllmhdims.pop(1)
464        cllmhvdims.pop(1)
465
466        cllmh = fdin.module_fordiagnostics.compute_cllmh4d2(cldfra[:], pres[:])
467       
468    else:
469        dx = cldfra.shape[2]
470        dy = cldfra.shape[1]
471        dz = cldfra.shape[0]
472        cllmhdims.pop(0)
473        cllmhvdims.pop(0)
474
475        cllmh = fdin.module_fordiagnostics.compute_cllmh3d1(cldfra[:], pres[:])
476
477    return cllmh, cllmhdims, cllmhvdims
478
479def compute_cllmh(cldfra, pres, dimns, dimvns):
480    """ Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ
481    compute_clt(cldfra, pres, dimns, dimvns)
482      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
483      [pres] = pressure field
484      [dimns]= list of the name of the dimensions of [cldfra]
485      [dimvns]= list of the name of the variables with the values of the
486        dimensions of [cldfra]
487    """
488    fname = 'compute_cllmh'
489
490    cllmhdims = dimns[:]
491    cllmhvdims = dimvns[:]
492
493    if len(cldfra.shape) == 4:
494        dx = cldfra.shape[3]
495        dy = cldfra.shape[2]
496        dz = cldfra.shape[1]
497        dt = cldfra.shape[0]
498        cllmhdims.pop(1)
499        cllmhvdims.pop(1)
500
501        cllmh = np.ones(tuple([3, dt, dy, dx]), dtype=np.float)
502
503        for it in range(dt):
504            for ix in range(dx):
505                for iy in range(dy):
506                    gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
507                    cllmh[:,it,iy,ix] = var_cllmh(cldfra[it,:,iy,ix], pres[it,:,iy,ix])
508       
509    else:
510        dx = cldfra.shape[2]
511        dy = cldfra.shape[1]
512        dz = cldfra.shape[0]
513        cllmhdims.pop(0)
514        cllmhvdims.pop(0)
515
516        cllmh = np.ones(tuple([3, dy, dx]), dtype=np.float)
517
518        for ix in range(dx):
519            for iy in range(dy):
520                gen.percendone(ix*dy + iy,dx*dy, 5, 'diagnosted')
521                cllmh[:,iy,ix] = var_cllmh(cldfra[:,iy,ix], pres[:,iy,ix])
522
523    return cllmh, cllmhdims, cllmhvdims
524
525def compute_clivi(dens, qtot, dimns, dimvns):
526    """ Function to compute cloud-ice water path (clivi)
527      [dens] = density [in kgkg-1] (assuming [t],z,y,x)
528      [qtot] = added mixing ratio of all cloud-ice species in [kgkg-1] (assuming [t],z,y,x)
529      [dimns]= list of the name of the dimensions of [q]
530      [dimvns]= list of the name of the variables with the values of the
531        dimensions of [q]
532    """
533    fname = 'compute_clivi'
534
535    clividims = dimns[:]
536    clivivdims = dimvns[:]
537
538    if len(qtot.shape) == 4:
539        clividims.pop(1)
540        clivivdims.pop(1)
541    else:
542        clividims.pop(0)
543        clivivdims.pop(0)
544
545    data1 = dens*qtot
546    clivi = np.sum(data1, axis=1)
547
548    return clivi, clividims, clivivdims
549
550
551def compute_clwvl(dens, qtot, dimns, dimvns):
552    """ Function to compute condensed water path (clwvl)
553      [dens] = density [in kgkg-1] (assuming [t],z,y,x)
554      [qtot] = added mixing ratio of all cloud-water species in [kgkg-1] (assuming [t],z,y,x)
555      [dimns]= list of the name of the dimensions of [q]
556      [dimvns]= list of the name of the variables with the values of the
557        dimensions of [q]
558    """
559    fname = 'compute_clwvl'
560
561    clwvldims = dimns[:]
562    clwvlvdims = dimvns[:]
563
564    if len(qtot.shape) == 4:
565        clwvldims.pop(1)
566        clwvlvdims.pop(1)
567    else:
568        clwvldims.pop(0)
569        clwvlvdims.pop(0)
570
571    data1 = dens*qtot
572    clwvl = np.sum(data1, axis=1)
573
574    return clwvl, clwvldims, clwvlvdims
575
576def var_virtualTemp (temp,rmix):
577    """ This function returns virtual temperature in K,
578      temp: temperature [K]
579      rmix: mixing ratio in [kgkg-1]
580    """
581
582    fname = 'var_virtualTemp'
583
584    virtual=temp*(0.622+rmix)/(0.622*(1.+rmix))
585
586    return virtual
587
588
589def var_mslp(pres, psfc, ter, tk, qv):
590    """ Function to compute mslp on a 1D column
591    """
592
593    fname = 'var_mslp'
594
595    N = 1.0
596    expon=287.04*.0065/9.81
597    pref = 40000.
598
599# First find where about 400 hPa is located
600    dz=len(pres) 
601
602    kref = -1
603    pinc = pres[0] - pres[dz-1]
604
605    if pinc < 0.:
606        for iz in range(1,dz):
607            if pres[iz-1] >= pref and pres[iz] < pref: 
608                kref = iz
609                break
610    else:
611        for iz in range(dz-1):
612            if pres[iz] >= pref and pres[iz+1] < pref: 
613                kref = iz
614                break
615
616    if kref == -1:
617        print errormsg
618        print '  ' + fname + ': no reference pressure:',pref,'found!!'
619        print '    values:',pres[:]
620        quit(-1)
621
622    mslp = 0.
623
624# We are below both the ground and the lowest data level.
625
626# First, find the model level that is closest to a "target" pressure
627# level, where the "target" pressure is delta-p less that the local
628# value of a horizontally smoothed surface pressure field.  We use
629# delta-p = 150 hPa here. A standard lapse rate temperature profile
630# passing through the temperature at this model level will be used
631# to define the temperature profile below ground.  This is similar
632# to the Benjamin and Miller (1990) method, using 
633# 700 hPa everywhere for the "target" pressure.
634
635# ptarget = psfc - 15000.
636    ptarget = 70000.
637    dpmin=1.e4
638    kupper = 0
639    if pinc > 0.:
640        for iz in range(dz-1,0,-1):
641            kupper = iz
642            dp=np.abs( pres[iz] - ptarget )
643            if dp < dpmin: exit
644            dpmin = np.min([dpmin, dp])
645    else:
646        for iz in range(dz):
647            kupper = iz
648            dp=np.abs( pres[iz] - ptarget )
649            if dp < dpmin: exit
650            dpmin = np.min([dpmin, dp])
651
652    pbot=np.max([pres[0], psfc])
653#    zbot=0.
654
655#    tbotextrap=tk(i,j,kupper,itt)*(pbot/pres_field(i,j,kupper,itt))**expon
656#    tvbotextrap=virtual(tbotextrap,qv(i,j,1,itt))
657
658#    data_out(i,j,itt,1) = (zbot+tvbotextrap/.0065*(1.-(interp_levels(1)/pbot)**expon))
659    tbotextrap = tk[kupper]*(psfc/ptarget)**expon
660    tvbotextrap = var_virtualTemp(tbotextrap, qv[kupper])
661    mslp = psfc*( (tvbotextrap+0.0065*ter)/tvbotextrap)**(1./expon)
662
663    return mslp
664
665def compute_mslp(pressure, psurface, terrain, temperature, qvapor, dimns, dimvns):
666    """ Function to compute mslp: mean sea level pressure following p_interp.F90 from WRF
667    var_mslp(pres, ter, tk, qv, dimns, dimvns)
668      [pressure]= pressure field [Pa] (assuming [[t],z,y,x])
669      [psurface]= surface pressure field [Pa]
670      [terrain]= topography [m]
671      [temperature]= temperature [K]
672      [qvapor]= water vapour mixing ratio [kgkg-1]
673      [dimns]= list of the name of the dimensions of [cldfra]
674      [dimvns]= list of the name of the variables with the values of the
675        dimensions of [pres]
676    """
677
678    fname = 'compute_mslp'
679
680    mslpdims = list(dimns[:])
681    mslpvdims = list(dimvns[:])
682
683    if len(pressure.shape) == 4:
684        mslpdims.pop(1)
685        mslpvdims.pop(1)
686    else:
687        mslpdims.pop(0)
688        mslpvdims.pop(0)
689
690    if len(pressure.shape) == 4:
691        dx = pressure.shape[3]
692        dy = pressure.shape[2]
693        dz = pressure.shape[1]
694        dt = pressure.shape[0]
695
696        mslpv = np.zeros(tuple([dt, dy, dx]), dtype=np.float)
697
698# Terrain... to 2D !
699        terval = np.zeros(tuple([dy, dx]), dtype=np.float)
700        if len(terrain.shape) == 3:
701            terval = terrain[0,:,:]
702        else:
703            terval = terrain
704
705        for ix in range(dx):
706            for iy in range(dy):
707                if terval[iy,ix] > 0.:
708                    for it in range(dt):
709                        mslpv[it,iy,ix] = var_mslp(pressure[it,:,iy,ix],             \
710                          psurface[it,iy,ix], terval[iy,ix], temperature[it,:,iy,ix],\
711                          qvapor[it,:,iy,ix])
712
713                        gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
714                else:
715                    mslpv[:,iy,ix] = psurface[:,iy,ix]
716
717    else:
718        dx = pressure.shape[2]
719        dy = pressure.shape[1]
720        dz = pressure.shape[0]
721
722        mslpv = np.zeros(tuple([dy, dx]), dtype=np.float)
723
724# Terrain... to 2D !
725        terval = np.zeros(tuple([dy, dx]), dtype=np.float)
726        if len(terrain.shape) == 3:
727            terval = terrain[0,:,:]
728        else:
729            terval = terrain
730
731        for ix in range(dx):
732            for iy in range(dy):
733                gen.percendone(ix*dy + iy,dx*dy, 5, 'diagnosted')
734                if terval[iy,ix] > 0.:
735                    mslpv[iy,ix] = var_mslp(pressure[:,iy,ix], psurface[iy,ix],          \
736                      terval[iy,ix], temperature[:,iy,ix], qvapor[:,iy,ix])
737                else:
738                    mslpv[iy,ix] = psfc[iy,ix]
739
740    return mslpv, mslpdims, mslpvdims
741
742def compute_OMEGAw(omega, p, t, dimns, dimvns):
743    """ Function to transform OMEGA [Pas-1] to velocities [ms-1]
744      tacking: https://www.ncl.ucar.edu/Document/Functions/Contributed/omega_to_w.shtml
745      [omega] = vertical velocity [in ms-1] (assuming [t],z,y,x)
746      [p] = pressure in [Pa] (assuming [t],z,y,x)
747      [t] = temperature in [K] (assuming [t],z,y,x)
748      [dimns]= list of the name of the dimensions of [q]
749      [dimvns]= list of the name of the variables with the values of the
750        dimensions of [q]
751    """
752    fname = 'compute_OMEGAw'
753
754    rgas = 287.058            # J/(kg-K) => m2/(s2 K)
755    g    = 9.80665            # m/s2
756
757    wdims = dimns[:]
758    wvdims = dimvns[:]
759
760    rho  = p/(rgas*t)         # density => kg/m3
761    w    = -omega/(rho*g)     
762
763    return w, wdims, wvdims
764
765def compute_prw(dens, q, dimns, dimvns):
766    """ Function to compute water vapour path (prw)
767      [dens] = density [in kgkg-1] (assuming [t],z,y,x)
768      [q] = mixing ratio in [kgkg-1] (assuming [t],z,y,x)
769      [dimns]= list of the name of the dimensions of [q]
770      [dimvns]= list of the name of the variables with the values of the
771        dimensions of [q]
772    """
773    fname = 'compute_prw'
774
775    prwdims = dimns[:]
776    prwvdims = dimvns[:]
777
778    if len(q.shape) == 4:
779        prwdims.pop(1)
780        prwvdims.pop(1)
781    else:
782        prwdims.pop(0)
783        prwvdims.pop(0)
784
785    data1 = dens*q
786    prw = np.sum(data1, axis=1)
787
788    return prw, prwdims, prwvdims
789
790def compute_rh(p, t, q, dimns, dimvns):
791    """ Function to compute relative humidity following 'Tetens' equation (T,P) ...'
792      [t]= temperature (assuming [[t],z,y,x] in [K])
793      [p] = pressure field (assuming in [hPa])
794      [q] = mixing ratio in [kgkg-1]
795      [dimns]= list of the name of the dimensions of [t]
796      [dimvns]= list of the name of the variables with the values of the
797        dimensions of [t]
798    """
799    fname = 'compute_rh'
800
801    rhdims = dimns[:]
802    rhvdims = dimvns[:]
803
804    data1 = 10.*0.6112*np.exp(17.67*(t-273.16)/(t-29.65))
805    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
806
807    rh = q/data2
808
809    return rh, rhdims, rhvdims
810
811def compute_td(p, temp, qv, dimns, dimvns):
812    """ Function to compute the dew point temperature
813      [p]= pressure [Pa]
814      [temp]= temperature [C]
815      [qv]= mixing ratio [kgkg-1]
816      [dimns]= list of the name of the dimensions of [p]
817      [dimvns]= list of the name of the variables with the values of the
818        dimensions of [p]
819    """
820    fname = 'compute_td'
821
822#    print '    ' + fname + ': computing dew-point temperature from TS as t and Tetens...'
823# tacking from: http://en.wikipedia.org/wiki/Dew_point
824    tk = temp
825    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
826    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
827
828    rh = qv/data2
829               
830    pa = rh * data1
831    td = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
832
833    tddims = dimns[:]
834    tdvdims = dimvns[:]
835
836    return td, tddims, tdvdims
837
838def var_WRFtime(timewrfv, refdate='19491201000000', tunitsval='minutes'):
839    """ Function to copmute CFtimes from WRFtime variable
840      refdate= [YYYYMMDDMIHHSS] format of reference date
841      tunitsval= CF time units
842      timewrfv= matrix string values of WRF 'Times' variable
843    """
844    fname = 'var_WRFtime'
845
846    yrref=refdate[0:4]
847    monref=refdate[4:6]
848    dayref=refdate[6:8]
849    horref=refdate[8:10]
850    minref=refdate[10:12]
851    secref=refdate[12:14]
852
853    refdateS = yrref + '-' + monref + '-' + dayref + ' ' + horref + ':' + minref +   \
854      ':' + secref
855
856    dt = timewrfv.shape[0]
857    WRFtime = np.zeros((dt), dtype=np.float)
858
859    for it in range(dt):
860        wrfdates = gen.datetimeStr_conversion(timewrfv[it,:],'WRFdatetime', 'matYmdHMS')
861        WRFtime[it] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
862
863    tunits = tunitsval + ' since ' + refdateS
864
865    return WRFtime, tunits
866
867def turbulence_var(varv, dimvn, dimn):
868    """ Function to compute the Taylor's decomposition turbulence term from a a given variable
869      x*=<x^2>_t-(<X>_t)^2
870    turbulence_var(varv,dimn)
871      varv= values of the variable
872      dimvn= names of the dimension of the variable
873      dimn= names of the dimensions (as a dictionary with 'X', 'Y', 'Z', 'T')
874    >>> turbulence_var(np.arange((27)).reshape(3,3,3),['time','y','x'],{'T':'time', 'Y':'y', 'X':'x'})
875    [[ 54.  54.  54.]
876     [ 54.  54.  54.]
877     [ 54.  54.  54.]]
878    """
879    fname = 'turbulence_varv'
880
881    timedimid = dimvn.index(dimn['T'])
882
883    varv2 = varv*varv
884
885    vartmean = np.mean(varv, axis=timedimid)
886    var2tmean = np.mean(varv2, axis=timedimid)
887
888    varvturb = var2tmean - (vartmean*vartmean)
889
890    return varvturb
891
892def compute_turbulence(v, dimns, dimvns):
893    """ Function to compute the rubulence term of the Taylor's decomposition ...'
894      x*=<x^2>_t-(<X>_t)^2
895      [v]= variable (assuming [[t],z,y,x])
896      [dimns]= list of the name of the dimensions of [v]
897      [dimvns]= list of the name of the variables with the values of the
898        dimensions of [v]
899    """
900    fname = 'compute_turbulence'
901
902    turbdims = dimns[:]
903    turbvdims = dimvns[:]
904
905    turbdims.pop(0)
906    turbvdims.pop(0)
907
908    v2 = v*v
909
910    vartmean = np.mean(v, axis=0)
911    var2tmean = np.mean(v2, axis=0)
912
913    turb = var2tmean - (vartmean*vartmean)
914
915    return turb, turbdims, turbvdims
916
917def compute_wds(u, v, dimns, dimvns):
918    """ Function to compute the wind direction
919      [u]= W-E wind direction [ms-1, knot, ...]
920      [v]= N-S wind direction [ms-1, knot, ...]
921      [dimns]= list of the name of the dimensions of [u]
922      [dimvns]= list of the name of the variables with the values of the
923        dimensions of [u]
924    """
925    fname = 'compute_wds'
926
927#    print '    ' + fname + ': computing wind direction as ATAN2(v,u) ...'
928    theta = np.arctan2(v,u)
929    theta = np.where(theta < 0., theta + 2.*np.pi, theta)
930
931    wds = 360.*theta/(2.*np.pi)
932
933    wdsdims = dimns[:]
934    wdsvdims = dimvns[:]
935
936    return wds, wdsdims, wdsvdims
937
938def compute_wss(u, v, dimns, dimvns):
939    """ Function to compute the wind speed
940      [u]= W-E wind direction [ms-1, knot, ...]
941      [v]= N-S wind direction [ms-1, knot, ...]
942      [dimns]= list of the name of the dimensions of [u]
943      [dimvns]= list of the name of the variables with the values of the
944        dimensions of [u]
945    """
946    fname = 'compute_wss'
947
948#    print '    ' + fname + ': computing wind speed as SQRT(v**2 + u**2) ...'
949    wss = np.sqrt(u*u + v*v)
950
951    wssdims = dimns[:]
952    wssvdims = dimvns[:]
953
954    return wss, wssdims, wssvdims
955
956def timeunits_seconds(dtu):
957    """ Function to transform a time units to seconds
958    timeunits_seconds(timeuv)
959      [dtu]= time units value to transform in seconds
960    """
961    fname='timunits_seconds'
962
963    if dtu == 'years':
964        times = 365.*24.*3600.
965    elif dtu == 'weeks':
966        times = 7.*24.*3600.
967    elif dtu == 'days':
968        times = 24.*3600.
969    elif dtu == 'hours':
970        times = 3600.
971    elif dtu == 'minutes':
972        times = 60.
973    elif dtu == 'seconds':
974        times = 1.
975    elif dtu == 'miliseconds':
976        times = 1./1000.
977    else:
978        print errormsg
979        print '  ' + fname  + ": time units '" + dtu + "' not ready !!"
980        quit(-1)
981
982    return times
983
984def compute_WRFua(u, v, sina, cosa, dimns, dimvns):
985    """ Function to compute geographical rotated WRF 3D winds
986      u= orginal WRF x-wind
987      v= orginal WRF y-wind
988      sina= original WRF local sinus of map rotation
989      cosa= original WRF local cosinus of map rotation
990      formula:
991        ua = u*cosa-va*sina
992        va = u*sina+va*cosa
993    """
994    fname = 'compute_WRFua'
995
996    var0 = u
997    var1 = v
998    var2 = sina
999    var3 = cosa
1000
1001    # un-staggering variables
1002    unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1003    ua = np.zeros(tuple(unstgdims), dtype=np.float)
1004    unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1005    unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1006    unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1007    unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1008
1009    for iz in range(var0.shape[1]):
1010        ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
1011
1012    dnamesvar = ['Time','bottom_top','south_north','west_east']
1013    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1014
1015    return ua, dnamesvar, dvnamesvar
1016
1017def compute_WRFva(u, v, sina, cosa, dimns, dimvns):
1018    """ Function to compute geographical rotated WRF 3D winds
1019      u= orginal WRF x-wind
1020      v= orginal WRF y-wind
1021      sina= original WRF local sinus of map rotation
1022      cosa= original WRF local cosinus of map rotation
1023      formula:
1024        ua = u*cosa-va*sina
1025        va = u*sina+va*cosa
1026    """
1027    fname = 'compute_WRFva'
1028
1029    var0 = u
1030    var1 = v
1031    var2 = sina
1032    var3 = cosa
1033
1034    # un-staggering variables
1035    unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1036    va = np.zeros(tuple(unstgdims), dtype=np.float)
1037    unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1038    unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1039    unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1040    unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1041
1042    for iz in range(var0.shape[1]):
1043        va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
1044
1045    dnamesvar = ['Time','bottom_top','south_north','west_east']
1046    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1047
1048    return va, dnamesvar, dvnamesvar
1049
1050def compute_WRFuava(u, v, sina, cosa, dimns, dimvns):
1051    """ Function to compute geographical rotated WRF 3D winds
1052      u= orginal WRF x-wind
1053      v= orginal WRF y-wind
1054      sina= original WRF local sinus of map rotation
1055      cosa= original WRF local cosinus of map rotation
1056      formula:
1057        ua = u*cosa-va*sina
1058        va = u*sina+va*cosa
1059    """
1060    fname = 'compute_WRFuava'
1061
1062    var0 = u
1063    var1 = v
1064    var2 = sina
1065    var3 = cosa
1066
1067    # un-staggering variables
1068    unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1069    ua = np.zeros(tuple(unstgdims), dtype=np.float)
1070    va = np.zeros(tuple(unstgdims), dtype=np.float)
1071    unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1072    unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1073    unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1074    unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1075
1076    for iz in range(var0.shape[1]):
1077        ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
1078        va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
1079
1080    dnamesvar = ['Time','bottom_top','south_north','west_east']
1081    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1082
1083    return ua, va, dnamesvar, dvnamesvar
1084
1085def compute_WRFuas(u10, v10, sina, cosa, dimns, dimvns):
1086    """ Function to compute geographical rotated WRF 2-meter x-wind
1087      u10= orginal WRF 10m x-wind
1088      v10= orginal WRF 10m y-wind
1089      sina= original WRF local sinus of map rotation
1090      cosa= original WRF local cosinus of map rotation
1091      formula:
1092        uas = u10*cosa-va10*sina
1093        vas = u10*sina+va10*cosa
1094    """
1095    fname = 'compute_WRFuas'
1096
1097    var0 = u10
1098    var1 = v10
1099    var2 = sina
1100    var3 = cosa
1101
1102    uas = np.zeros(var0.shape, dtype=np.float)
1103    vas = np.zeros(var0.shape, dtype=np.float)
1104
1105    uas = var0*var3 - var1*var2
1106
1107    dnamesvar = ['Time','south_north','west_east']
1108    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1109
1110    return uas, dnamesvar, dvnamesvar
1111
1112def compute_WRFvas(u10, v10, sina, cosa, dimns, dimvns):
1113    """ Function to compute geographical rotated WRF 2-meter y-wind
1114      u10= orginal WRF 10m x-wind
1115      v10= orginal WRF 10m y-wind
1116      sina= original WRF local sinus of map rotation
1117      cosa= original WRF local cosinus of map rotation
1118      formula:
1119        uas = u10*cosa-va10*sina
1120        vas = u10*sina+va10*cosa
1121    """
1122    fname = 'compute_WRFvas'
1123
1124    var0 = u10
1125    var1 = v10
1126    var2 = sina
1127    var3 = cosa
1128
1129    uas = np.zeros(var0.shape, dtype=np.float)
1130    vas = np.zeros(var0.shape, dtype=np.float)
1131
1132    vas = var0*var2 + var1*var3
1133
1134    dnamesvar = ['Time','south_north','west_east']
1135    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1136
1137    return vas, dnamesvar, dvnamesvar
1138
1139def compute_WRFuasvas(u10, v10, sina, cosa, dimns, dimvns):
1140    """ Function to compute geographical rotated WRF 2-meter winds
1141      u10= orginal WRF 10m x-wind
1142      v10= orginal WRF 10m y-wind
1143      sina= original WRF local sinus of map rotation
1144      cosa= original WRF local cosinus of map rotation
1145      formula:
1146        uas = u10*cosa-va10*sina
1147        vas = u10*sina+va10*cosa
1148    """
1149    fname = 'compute_WRFuasvas'
1150
1151    var0 = u10
1152    var1 = v10
1153    var2 = sina
1154    var3 = cosa
1155
1156    uas = np.zeros(var0.shape, dtype=np.float)
1157    vas = np.zeros(var0.shape, dtype=np.float)
1158
1159    uas = var0*var3 - var1*var2
1160    vas = var0*var2 + var1*var3
1161
1162    dnamesvar = ['Time','south_north','west_east']
1163    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1164
1165    return uas, vas, dnamesvar, dvnamesvar
1166
1167def compute_WRFta(t, p, dimns, dimvns):
1168    """ Function to compute WRF air temperature
1169      t= orginal WRF temperature
1170      p= original WRF pressure (P + PB)
1171      formula:
1172          temp = theta*(p/p0)**(R/Cp)
1173
1174    """
1175    fname = 'compute_WRFta'
1176
1177    ta = (t+300.)*(p/fdef.module_definitions.p0ref)**(fdef.module_definitions.rcp)
1178
1179    dnamesvar = ['Time','bottom_top','south_north','west_east']
1180    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1181
1182    return ta, dnamesvar, dvnamesvar
1183
1184def compute_WRFtd(t, p, qv, dimns, dimvns):
1185    """ Function to compute WRF dew-point air temperature
1186      t= orginal WRF temperature
1187      p= original WRF pressure (P + PB)
1188      formula:
1189          temp = theta*(p/p0)**(R/Cp)
1190
1191    """
1192    fname = 'compute_WRFtd'
1193
1194    tk = (t+300.)*(p/fdef.module_definitions.p0ref)**(fdef.module_definitions.rcp)
1195    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
1196    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
1197
1198    rh = qv/data2
1199               
1200    pa = rh * data1
1201    td = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
1202
1203    dnamesvar = ['Time','bottom_top','south_north','west_east']
1204    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1205
1206    return td, dnamesvar, dvnamesvar
1207
1208def compute_WRFwd(u, v, sina, cosa, dimns, dimvns):
1209    """ Function to compute the wind direction
1210      u= W-E wind direction [ms-1]
1211      v= N-S wind direction [ms-1]
1212      sina= original WRF local sinus of map rotation
1213      cosa= original WRF local cosinus of map rotation
1214    """
1215    fname = 'compute_WRFwd'
1216    var0 = u
1217    var1 = v
1218    var2 = sina
1219    var3 = cosa
1220
1221    # un-staggering variables
1222    unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1223    ua = np.zeros(tuple(unstgdims), dtype=np.float)
1224    va = np.zeros(tuple(unstgdims), dtype=np.float)
1225    unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1226    unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1227    unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1228    unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1229
1230    for iz in range(var0.shape[1]):
1231        ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
1232        va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
1233
1234    theta = np.arctan2(va,ua)
1235    theta = np.where(theta < 0., theta + 2.*np.pi, theta)
1236
1237    wd = 360.*theta/(2.*np.pi)
1238
1239    dnamesvar = ['Time','bottom_top','south_north','west_east']
1240    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1241
1242    return wd
1243
1244def var_td(t, p, qv):
1245    """ Function to compute dew-point air temperature from temperature and pressure values
1246      t= temperature [K]
1247      p= pressure (Pa)
1248      formula:
1249          temp = theta*(p/p0)**(R/Cp)
1250
1251    """
1252    fname = 'compute_td'
1253
1254    tk = (t)*(p/fdef.module_definitions.p0ref)**(fdef.module_definitions.rcp)
1255    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
1256    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
1257
1258    rh = qv/data2
1259               
1260    pa = rh * data1
1261    td = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
1262
1263    return td
1264
1265def var_wd(u, v):
1266    """ Function to compute the wind direction
1267      [u]= W-E wind direction [ms-1, knot, ...]
1268      [v]= N-S wind direction [ms-1, knot, ...]
1269    """
1270    fname = 'var_wd'
1271
1272    theta = np.arctan2(v,u)
1273    theta = np.where(theta < 0., theta + 2.*np.pi, theta)
1274
1275    wd = 360.*theta/(2.*np.pi)
1276
1277    return wd
1278
1279def var_ws(u, v):
1280    """ Function to compute the wind speed
1281      [u]= W-E wind direction [ms-1, knot, ...]
1282      [v]= N-S wind direction [ms-1, knot, ...]
1283    """
1284    fname = 'var_ws'
1285
1286    ws = np.sqrt(u*u + v*v)
1287
1288    return ws
1289
1290class C_diagnostic(object):
1291    """ Class to compute generic variables
1292      Cdiag: name of the diagnostic to compute
1293      ncobj: netcdf object with data
1294      sfcvars: dictionary with CF equivalencies of surface variables inside file
1295      vars3D: dictionary with CF equivalencies of 3D variables inside file
1296      dictdims: dictionary with CF equivalencies of dimensions inside file
1297        self.values = Values of the diagnostic
1298        self.dims = Dimensions of the diagnostic
1299        self.units = units of the diagnostic
1300        self.incvars = list of variables from the input netCDF object
1301    """ 
1302    def __init__(self, Cdiag, ncobj, sfcvars, vars3D, dictdims):
1303        fname = 'C_diagnostic'
1304        self.values = None
1305        self.dims = None
1306        self.incvars = ncobj.variables
1307        self.units = None
1308
1309        if Cdiag == 'td':
1310            """ Computing dew-point temperature
1311            """
1312            vn = 'td'
1313            CF3Dvars = ['ta', 'plev', 'hus']
1314            for v3D in CF3Dvars:
1315                if not vars3D.has_key(v3D):
1316                    print gen.errormsg
1317                    print '  ' + fname + ": missing variable '" + v3D +              \
1318                      "' attribution to compute '" + vn + "' !!"
1319                    print '  Equivalence of 3D variables provided _______'
1320                    gen.printing_dictionary(vars3D)
1321                    quit(-1)
1322                if not self.incvars.has_key(vars3D[v3D]):
1323                    print gen.errormsg
1324                    print '  ' + fname + ": missing variable '" + vars3D[v3D] +      \
1325                      "' in input file to compute '" + vn + "' !!"
1326                    print '  available variables:', self.incvars.keys()
1327                    print '  looking for variables _______' 
1328                    gen.printing_dictionary(vars3D)
1329                    quit(-1)
1330
1331            ta = ncobj.variables[vars3D['ta']][:]
1332            p = ncobj.variables[vars3D['plev']][:]
1333            hur = ncobj.variables[vars3D['hus']][:]
1334
1335            if len(ta.shape) != len(p.shape):
1336                p = fill_Narray(p, ta*0., filldim=[0,2,3])
1337
1338            self.values = var_td(ta, p, hur)
1339            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1340              dictdims['lon']]
1341            self.units = 'K'
1342
1343        elif Cdiag == 'wd':
1344            """ Computing wind direction
1345            """
1346            vn = 'wd'
1347            CF3Dvars = ['ua', 'va']
1348            for v3D in CF3Dvars:
1349                if not vars3D.has_key(v3D):
1350                    print gen.errormsg
1351                    print '  ' + fname + ": missing variable '" + v3D +              \
1352                      "self.' attribution to compute '" + vn + "' !!"
1353                    print '  Equivalence of 3D variables provided _______'
1354                    gen.printing_dictionary(vars3D)
1355                    quit(-1)
1356                if not self.incvars.has_key(vars3D[v3D]):
1357                    print gen.errormsg
1358                    print '  ' + fname + ": missing variable '" + vars3D[v3D] +      \
1359                      "' in input file to compute '" + vn + "' !!"
1360                    print '  available variables:', self.incvars.keys()
1361                    print '  looking for variables _______' 
1362                    gen.printing_dictionary(vars3D)
1363                    quit(-1)
1364   
1365            ua = ncobj.variables[vars3D['ua']][:]
1366            va = ncobj.variables[vars3D['va']][:]
1367   
1368            self.values = var_wd(ua, va)
1369            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1370              dictdims['lon']]
1371            self.units = 'degree'
1372
1373        elif Cdiag == 'ws':
1374            """ Computing wind speed
1375            """
1376            vn = 'ws'
1377            CF3Dvars = ['ua', 'va']
1378            for v3D in CF3Dvars:
1379                if not vars3D.has_key(v3D):
1380                    print gen.errormsg
1381                    print '  ' + fname + ": missing variable '" + v3D +              \
1382                      "' attribution to compute '" + vn + "' !!"
1383                    print '  Equivalence of 3D variables provided _______'
1384                    gen.printing_dictionary(vars3D)
1385                    quit(-1)
1386                if not self.incvars.has_key(vars3D[v3D]):
1387                    print gen.errormsg
1388                    print '  ' + fname + ": missing variable '" + vars3D[v3D] +      \
1389                      "' in input file to compute '" + vn + "' !!"
1390                    print '  available variables:', self.incvars.keys()
1391                    print '  looking for variables _______' 
1392                    gen.printing_dictionary(vars3D)
1393                    quit(-1)
1394
1395            ua = ncobj.variables[vars3D['ua']][:]
1396            va = ncobj.variables[vars3D['va']][:]
1397
1398            self.values = var_ws(ua, va)
1399            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1400              dictdims['lon']]
1401            self.units = ncobj.variables[vars3D['ua']].units
1402
1403        else:
1404            print gen.errormsg
1405            print '  ' + fname + ": variable '" + Wdiag + "' not ready !!"
1406            print '  available ones:', Cavailablediags
1407            quit(-1)
1408
1409class W_diagnostic(object):
1410    """ Class to compute WRF diagnostics variables
1411      Wdiag: name of the diagnostic to compute
1412      ncobj: netcdf object with data
1413      sfcvars: dictionary with CF equivalencies of surface variables inside file
1414      vars3D: dictionary with CF equivalencies of 3D variables inside file
1415      indims: list of dimensions inside file
1416      invardims: list of dimension-variables inside file
1417      dictdims: dictionary with CF equivalencies of dimensions inside file
1418        self.values = Values of the diagnostic
1419        self.dims = Dimensions of the diagnostic
1420        self.units = units of the diagnostic
1421        self.incvars = list of variables from the input netCDF object
1422    """   
1423    def __init__(self, Wdiag, ncobj, sfcvars, vars3D, indims, invardims, dictdims):
1424        fname = 'W_diagnostic'
1425
1426        self.values = None
1427        self.dims = None
1428        self.incvars = ncobj.variables
1429        self.units = None
1430
1431        if Wdiag == 'p':
1432            """ Computing air pressure
1433            """
1434            vn = 'p'
1435
1436            self.values = ncobj.variables['PB'][:] + ncobj.variables['P'][:]
1437            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1438              dictdims['lon']]
1439            self.units = ncobj.variables['PB'].units
1440
1441        elif Wdiag == 'ta':
1442            """ Computing air temperature
1443            """
1444            vn = 'ta'
1445            CF3Dvars = ['ta']
1446            for v3D in CF3Dvars:
1447                if not vars3D.has_key(v3D):
1448                    print gen.errormsg
1449                    print '  ' + fname + ": missing variable '" + v3D +              \
1450                      "' attribution to compute '" + vn + "' !!"
1451                    print '  Equivalence of 3D variables provided _______'
1452                    gen.printing_dictionary(vars3D)
1453                    quit(-1)
1454
1455            ta = ncobj.variables['T'][:]
1456            p = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
1457   
1458            vals, dims, vdims = compute_WRFta(ta, p, indims, invardims)
1459            self.values = vals
1460            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1461              dictdims['lon']]
1462            self.units = 'K'
1463
1464        elif Wdiag == 'td':
1465            """ Computing dew-point temperature
1466            """
1467            vn = 'td'
1468            CF3Dvars = ['ta', 'hus']
1469            for v3D in CF3Dvars:
1470                if not vars3D.has_key(v3D):
1471                    print gen.errormsg
1472                    print '  ' + fname + ": missing variable '" + v3D +              \
1473                      "' attribution to compute '" + vn + "' !!"
1474                    print '  Equivalence of 3D variables provided _______'
1475                    gen.printing_dictionary(vars3D)
1476                    quit(-1)
1477
1478            ta = ncobj.variables['T'][:]
1479            p = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
1480            hur = ncobj.variables['QVAPOR'][:]
1481   
1482            vals, dims, vdims = compute_WRFtd(ta, p, hur, indims, invardims)
1483            self.values = vals
1484            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1485              dictdims['lon']]
1486            self.units = 'K'
1487   
1488        elif Wdiag == 'ua':
1489            """ Computing x-wind
1490            """
1491            vn = 'ua'
1492            CF3Dvars = ['ua', 'va']
1493            for v3D in CF3Dvars:
1494                if not vars3D.has_key(v3D):
1495                    print gen.errormsg
1496                    print '  ' + fname + ": missing variable '" + v3D +              \
1497                      "' attribution to compute '" + vn + "' !!"
1498                    print '  Equivalence of 3D variables provided _______'
1499                    gen.printing_dictionary(vars3D)
1500                    quit(-1)
1501
1502            ua = ncobj.variables['U'][:]
1503            va = ncobj.variables['V'][:]
1504            sina = ncobj.variables['SINALPHA'][:]
1505            cosa = ncobj.variables['COSALPHA'][:]
1506   
1507            vals, dims, vdims = compute_WRFua(ua, va, sina, cosa, indims, invardims)
1508            self.values = vals
1509            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1510              dictdims['lon']]
1511            self.units = ncobj.variables['U'].units
1512   
1513        elif Wdiag == 'uas':
1514            """ Computing 10m x-wind
1515            """
1516            vn = 'uas'
1517            CFsfcvars = ['uas', 'vas']
1518            for vsf in CFsfcvars:
1519                if not sfcvars.has_key(vsf):
1520                    print gen.errormsg
1521                    print '  ' + fname + ": missing variable '" + vsf +              \
1522                      "' attribution to compute '" + vn + "' !!"
1523                    print '  Equivalence of sfc variables provided _______'
1524                    gen.printing_dictionary(sfcvars)
1525                    quit(-1)
1526   
1527            uas = ncobj.variables['U10'][:]
1528            vas = ncobj.variables['V10'][:]
1529            sina = ncobj.variables['SINALPHA'][:]
1530            cosa = ncobj.variables['COSALPHA'][:]
1531   
1532            vals,dims,vdims = compute_WRFuas(uas, vas, sina, cosa, indims, invardims)
1533            self.values = vals
1534            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1535              dictdims['lon']]
1536            self.units = ncobj.variables['U10'].units
1537   
1538        elif Wdiag == 'va':
1539            """ Computing y-wind
1540            """
1541            vn = 'ua'
1542            CF3Dvars = ['ua', 'va']
1543            for v3D in CF3Dvars:
1544                if not vars3D.has_key(v3D):
1545                    print gen.errormsg
1546                    print '  ' + fname + ": missing variable '" + v3D +              \
1547                      "' attribution to compute '" + vn + "' !!"
1548                    print '  Equivalence of 3D variables provided _______'
1549                    gen.printing_dictionary(vars3D)
1550                    quit(-1)
1551   
1552            ua = ncobj.variables['U'][:]
1553            va = ncobj.variables['V'][:]
1554            sina = ncobj.variables['SINALPHA'][:]
1555            cosa = ncobj.variables['COSALPHA'][:]
1556   
1557            vals, dims, vdims = compute_WRFva(ua, va, sina, cosa, indims, invardims)
1558            self.values = vals
1559            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1560              dictdims['lon']]
1561            self.units = ncobj.variables['U'].units
1562   
1563        elif Wdiag == 'vas':
1564            """ Computing 10m y-wind
1565            """
1566            vn = 'uas'
1567            CFsfcvars = ['uas', 'vas']
1568            for vsf in CFsfcvars:
1569                if not sfcvars.has_key(vsf):
1570                    print gen.errormsg
1571                    print '  ' + fname + ": missing variable '" + vsf +              \
1572                      "' attribution to compute '" + vn + "' !!"
1573                    print '  Equivalence of sfc variables provided _______'
1574                    gen.printing_dictionary(sfcvars)
1575                    quit(-1)
1576   
1577            uas = ncobj.variables['U10'][:]
1578            vas = ncobj.variables['V10'][:]
1579            sina = ncobj.variables['SINALPHA'][:]
1580            cosa = ncobj.variables['COSALPHA'][:]
1581   
1582            vals,dims,vdims = compute_WRFvas(uas, vas, sina, cosa, indims, invardims)
1583            self.values = vals
1584            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1585              dictdims['lon']]
1586            self.units = ncobj.variables['U10'].units
1587
1588        elif Wdiag == 'wd':
1589            """ Computing wind direction
1590            """
1591            vn = 'wd'
1592            CF3Dvars = ['ua', 'va']
1593            for v3D in CF3Dvars:
1594                if not vars3D.has_key(v3D):
1595                    print gen.errormsg
1596                    print '  ' + fname + ": missing variable '" + v3D +              \
1597                      "' attribution to compute '" + vn + "' !!"
1598                    print '  Equivalence of 3D variables provided _______'
1599                    gen.printing_dictionary(vars3D)
1600                    quit(-1)
1601
1602            ua = ncobj.variables['U10'][:]
1603            va = ncobj.variables['V10'][:]
1604            sina = ncobj.variables['SINALPHA'][:]
1605            cosa = ncobj.variables['COSALPHA'][:]
1606   
1607            vals, dims, vdims = compute_WRFwd(ua, va, sina, cosa, indims, invardims)
1608            self.values = vals
1609            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1610              dictdims['lon']]
1611            self.units = 'degree'
1612
1613        elif Wdiag == 'ws':
1614            """ Computing wind speed
1615            """
1616            vn = 'ws'
1617            CF3Dvars = ['ua', 'va']
1618            for v3D in CF3Dvars:
1619                if not vars3D.has_key(v3D):
1620                    print gen.errormsg
1621                    print '  ' + fname + ": missing variable '" + v3D +              \
1622                      "' attribution to compute '" + vn + "' !!"
1623                    print '  Equivalence of 3D variables provided _______'
1624                    gen.printing_dictionary(vars3D)
1625                    quit(-1)
1626   
1627            ua = ncobj.variables['U10'][:]
1628            va = ncobj.variables['V10'][:]
1629   
1630            self.values = var_ws(ua, va)
1631            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1632              dictdims['lon']]
1633            self.units = ncobj.variables['U10'].units
1634
1635        elif Wdiag == 'zg':
1636            """ Computing geopotential
1637            """
1638            vn = 'zg'
1639
1640            self.values = ncobj.variables['PHB'][:] + ncobj.variables['PH'][:]
1641            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1642              dictdims['lon']]
1643            self.units = ncobj.variables['PHB'].units
1644
1645        else:
1646            print gen.errormsg
1647            print '  ' + fname + ": variable '" + Wdiag + "' not ready !!"
1648            print '  available ones:', Wavailablediags
1649            quit(-1)
Note: See TracBrowser for help on using the repository browser.