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

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

Adding:

`psl_ecmwf': sea-level pressure computation following ECMWF

File size: 78.1 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_WRFhur: Function to compute WRF relative humidity following Teten's equation
22# compute_WRFta: Function to compute WRF air temperature
23# compute_WRFtd: Function to compute WRF dew-point air temperature
24# compute_WRFua: Function to compute geographical rotated WRF x-wind
25# compute_WRFva: Function to compute geographical rotated WRF y-wind
26# compute_WRFuava: Function to compute geographical rotated WRF 3D winds
27# compute_WRFuas: Function to compute geographical rotated WRF 2-meter x-wind
28# compute_WRFvas: Function to compute geographical rotated WRF 2-meter y-wind
29# compute_WRFuasvas: Fucntion to compute geographical rotated WRF 2-meter winds
30# derivate_centered: Function to compute the centered derivate of a given field
31# Forcompute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ via Fortran subroutine
32# Forcompute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ via a Fortran module
33# Forcompute_psl_ptarget: Function to compute the sea-level pressure following target_pressure value found in `p_interp.F'
34# Forcompute_zmla_gen: Function to compute the boundary layer height following a generic method with Fortran
35# Forcompute_zwind: Function to compute the wind at a given height following the power law method
36# Forcompute_zwind_log: Function to compute the wind at a given height following the logarithmic law method
37# Forcompute_zwindMO: Function to compute the wind at a given height following the Monin-Obukhov theory
38# W_diagnostic: Class to compute WRF diagnostics variables
39
40# Others just providing variable values
41# var_cllmh: Fcuntion to compute cllmh on a 1D column
42# var_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ using 1D vertical column values
43# var_hur: Function to compute relative humidity following 'August - Roche - Magnus' formula
44# var_hur_Uhus: Function to compute relative humidity following 'August - Roche - Magnus' formula using hus
45# var_mslp: Fcuntion to compute mean sea-level pressure
46# var_td: Function to compute dew-point air temperature from temperature and pressure values
47# var_td_Uhus: Function to compute dew-point air temperature from temperature and pressure values using hus
48# var_virtualTemp: This function returns virtual temperature in K,
49# var_WRFtime: Function to copmute CFtimes from WRFtime variable
50# var_wd: Function to compute the wind direction
51# var_wd: Function to compute the wind speed
52# rotational_z: z-component of the rotatinoal of horizontal vectorial field
53# turbulence_var: Function to compute the Taylor's decomposition turbulence term from a a given variable
54
55import numpy as np
56from netCDF4 import Dataset as NetCDFFile
57import os
58import re
59import nc_var_tools as ncvar
60import generic_tools as gen
61import datetime as dtime
62import module_ForDiag as fdin
63import module_ForDef as fdef
64
65main = 'diag_tools.py'
66errormsg = 'ERROR -- error -- ERROR -- error'
67warnmsg = 'WARNING -- warning -- WARNING -- warning'
68
69# Constants
70grav = fdef.module_definitions.grav
71
72# Available WRFiag
73Wavailablediags = ['hur', 'p', 'ta', 'td', 'ua', 'va', 'uas', 'vas', 'wd', 'ws', 'zg']
74
75# Available General diagnostics
76Cavailablediags = ['hur', 'hur_Uhus', 'td', 'td_Uhus', 'wd', 'ws']
77
78# Gneral information
79##
80def reduce_spaces(string):
81    """ Function to give words of a line of text removing any extra space
82    """
83    values = string.replace('\n','').split(' ')
84    vals = []
85    for val in values:
86         if len(val) > 0:
87             vals.append(val)
88
89    return vals
90
91def variable_combo(varn,combofile):
92    """ Function to provide variables combination from a given variable name
93      varn= name of the variable
94      combofile= ASCII file with the combination of variables
95        [varn] [combo]
96          [combo]: '@' separated list of variables to use to generate [varn]
97            [WRFdt] to get WRF time-step (from general attributes)
98    >>> variable_combo('WRFprls','/home/lluis/PY/diagnostics.inf')
99    deaccum@RAINNC@XTIME@prnc
100    """
101    fname = 'variable_combo'
102
103    if varn == 'h':
104        print fname + '_____________________________________________________________'
105        print variable_combo.__doc__
106        quit()
107
108    if not os.path.isfile(combofile):
109        print errormsg
110        print '  ' + fname + ": file with combinations '" + combofile +              \
111          "' does not exist!!"
112        quit(-1)
113
114    objf = open(combofile, 'r')
115
116    found = False
117    for line in objf:
118        linevals = reduce_spaces(line)
119        varnf = linevals[0]
120        combo = linevals[1].replace('\n','')
121        if varn == varnf: 
122            found = True
123            break
124
125    if not found:
126        print errormsg
127        print '  ' + fname + ": variable '" + varn + "' not found in '" + combofile +\
128          "' !!"
129        combo='ERROR'
130
131    objf.close()
132
133    return combo
134
135# Mathematical operators
136##
137def compute_accum(varv, dimns, dimvns):
138    """ Function to compute the accumulation of a variable
139    compute_accum(varv, dimnames, dimvns)
140      [varv]= values to accum (assuming [t,])
141      [dimns]= list of the name of the dimensions of the [varv]
142      [dimvns]= list of the name of the variables with the values of the
143        dimensions of [varv]
144    """
145    fname = 'compute_accum'
146
147    deacdims = dimns[:]
148    deacvdims = dimvns[:]
149
150    slicei = []
151    slicee = []
152
153    Ndims = len(varv.shape)
154    for iid in range(0,Ndims):
155        slicei.append(slice(0,varv.shape[iid]))
156        slicee.append(slice(0,varv.shape[iid]))
157
158    slicee[0] = np.arange(varv.shape[0])
159    slicei[0] = np.arange(varv.shape[0])
160    slicei[0][1:varv.shape[0]] = np.arange(varv.shape[0]-1)
161
162    vari = varv[tuple(slicei)]
163    vare = varv[tuple(slicee)]
164
165    ac = vari*0.
166    for it in range(1,varv.shape[0]):
167        ac[it,] = ac[it-1,] + vare[it,]
168
169    return ac, deacdims, deacvdims
170
171def compute_deaccum(varv, dimns, dimvns):
172    """ Function to compute the deaccumulation of a variable
173    compute_deaccum(varv, dimnames, dimvns)
174      [varv]= values to deaccum (assuming [t,])
175      [dimns]= list of the name of the dimensions of the [varv]
176      [dimvns]= list of the name of the variables with the values of the
177        dimensions of [varv]
178    """
179    fname = 'compute_deaccum'
180
181    deacdims = dimns[:]
182    deacvdims = dimvns[:]
183
184    slicei = []
185    slicee = []
186
187    Ndims = len(varv.shape)
188    for iid in range(0,Ndims):
189        slicei.append(slice(0,varv.shape[iid]))
190        slicee.append(slice(0,varv.shape[iid]))
191
192    slicee[0] = np.arange(varv.shape[0])
193    slicei[0] = np.arange(varv.shape[0])
194    slicei[0][1:varv.shape[0]] = np.arange(varv.shape[0]-1)
195
196    vari = varv[tuple(slicei)]
197    vare = varv[tuple(slicee)]
198
199    deac = vare - vari
200
201    return deac, deacdims, deacvdims
202
203def derivate_centered(var,dim,dimv):
204    """ Function to compute the centered derivate of a given field
205      centered derivate(n) = (var(n-1) + var(n+1))/(2*dn).
206    [var]= variable
207    [dim]= which dimension to compute the derivate
208    [dimv]= dimension values (can be of different dimension of [var])
209    >>> derivate_centered(np.arange(16).reshape(4,4)*1.,1,1.)
210    [[  0.   1.   2.   0.]
211     [  0.   5.   6.   0.]
212     [  0.   9.  10.   0.]
213     [  0.  13.  14.   0.]]
214    """
215
216    fname = 'derivate_centered'
217   
218    vark = var.dtype
219
220    if hasattr(dimv, "__len__"):
221# Assuming that the last dimensions of var [..., N, M] are the same of dimv [N, M]
222        if len(var.shape) != len(dimv.shape):
223            dimvals = np.zeros((var.shape), dtype=vark)
224            if len(var.shape) - len(dimv.shape) == 1:
225                for iz in range(var.shape[0]):
226                    dimvals[iz,] = dimv
227            elif len(var.shape) - len(dimv.shape) == 2:
228                for it in range(var.shape[0]):
229                    for iz in range(var.shape[1]):
230                        dimvals[it,iz,] = dimv
231            else:
232                print errormsg
233                print '  ' + fname + ': dimension difference between variable',      \
234                  var.shape,'and variable with dimension values',dimv.shape,         \
235                  ' not ready !!!'
236                quit(-1)
237        else:
238            dimvals = dimv
239    else:
240# dimension values are identical everywhere!
241# from: http://stackoverflow.com/questions/16807011/python-how-to-identify-if-a-variable-is-an-array-or-a-scalar   
242        dimvals = np.ones((var.shape), dtype=vark)*dimv
243
244    derivate = np.zeros((var.shape), dtype=vark)
245    if dim > len(var.shape) - 1:
246        print errormsg
247        print '  ' + fname + ': dimension',dim,' too big for given variable of ' +   \
248          'shape:', var.shape,'!!!'
249        quit(-1)
250
251    slicebef = []
252    sliceaft = []
253    sliceder = []
254
255    for id in range(len(var.shape)):
256        if id == dim:
257            slicebef.append(slice(0,var.shape[id]-2))
258            sliceaft.append(slice(2,var.shape[id]))
259            sliceder.append(slice(1,var.shape[id]-1))
260        else:
261            slicebef.append(slice(0,var.shape[id]))
262            sliceaft.append(slice(0,var.shape[id]))
263            sliceder.append(slice(0,var.shape[id]))
264
265    if hasattr(dimv, "__len__"):
266        derivate[tuple(sliceder)] = (var[tuple(slicebef)] + var[tuple(sliceaft)])/   \
267          ((dimvals[tuple(sliceaft)] - dimvals[tuple(slicebef)]))
268        print (dimvals[tuple(sliceaft)] - dimvals[tuple(slicebef)])
269    else:
270        derivate[tuple(sliceder)] = (var[tuple(slicebef)] + var[tuple(sliceaft)])/   \
271          (2.*dimv)
272
273#    print 'before________'
274#    print var[tuple(slicebef)]
275
276#    print 'after________'
277#    print var[tuple(sliceaft)]
278
279    return derivate
280
281def rotational_z(Vx,Vy,pos):
282    """ z-component of the rotatinoal of horizontal vectorial field
283    \/ x (Vx,Vy,Vz) = \/xVy - \/yVx
284    [Vx]= Variable component x
285    [Vy]=  Variable component y
286    [pos]= poisition of the grid points
287    >>> rotational_z(np.arange(16).reshape(4,4)*1., np.arange(16).reshape(4,4)*1., 1.)
288    [[  0.   1.   2.   0.]
289     [ -4.   0.   0.  -7.]
290     [ -8.   0.   0. -11.]
291     [  0.  13.  14.   0.]]
292    """
293
294    fname =  'rotational_z'
295
296    ndims = len(Vx.shape)
297    rot1 = derivate_centered(Vy,ndims-1,pos)
298    rot2 = derivate_centered(Vx,ndims-2,pos)
299
300    rot = rot1 - rot2
301
302    return rot
303
304# Diagnostics
305##
306def Forcompute_cape_afwa(ta, hur, pa, zg, hgt, parcelm, dimns, dimvns):
307    """ Function to compute the CAPE, CIN, ZLFC, PLFC, LI following WRF
308      'phys/module_diaf_afwa.F' methodology
309    Forcompute_cape_afwa(ta, hur, pa, hgt, zsfc, parcelm, dimns, dimvns)
310      [ta]= air-temperature values (assuming [[t],z,y,x]) [K]
311      [pa]= pressure values (assuming [[t],z,y,x]) [Pa]
312      [zg]= gopotential height (assuming [[t],z,y,x]) [gpm]
313      [hgt]= topographical height (assuming [m]
314      [parcelm]= method of air-parcel to use
315      [dimns]= list of the name of the dimensions of [pa]
316      [dimvns]= list of the name of the variables with the values of the
317        dimensions of [pa]
318    """
319    fname = 'Forcompute_cape_afwa'
320
321    psldims = dimns[:]
322    pslvdims = dimvns[:]
323
324    if len(pa.shape) == 4:
325        cape = np.zeros((pa.shape[0],pa.shape[2],pa.shape[3]), dtype=np.float)
326        cin = np.zeros((pa.shape[0],pa.shape[2],pa.shape[3]), dtype=np.float)
327        zlfc = np.zeros((pa.shape[0],pa.shape[2],pa.shape[3]), dtype=np.float)
328        plfc = np.zeros((pa.shape[0],pa.shape[2],pa.shape[3]), dtype=np.float)
329        li = np.zeros((pa.shape[0],pa.shape[2],pa.shape[3]), dtype=np.float)
330
331        dx = pa.shape[3]
332        dy = pa.shape[2]
333        dz = pa.shape[1]
334        dt = pa.shape[0]
335        psldims.pop(1)
336        pslvdims.pop(1)
337
338        pcape,pcin,pzlfc,pplfc,pli= fdin.module_fordiagnostics.compute_cape_afwa4d(  \
339          ta=ta[:].transpose(), hur=hur[:].transpose(), press=pa[:].transpose(),     \
340          zg=zg[:].transpose(), hgt=hgt.transpose(), parcelmethod=parcelm, d1=dx,    \
341          d2=dy, d3=dz, d4=dt)
342        cape = pcape.transpose()
343        cin = pcin.transpose()
344        zlfc = pzlfc.transpose()
345        plfc = pplfc.transpose()
346        li = pli.transpose()
347    else:
348        print errormsg
349        print '  ' + fname + ': rank', len(pa.shape), 'not ready !!'
350        print '  it only computes 4D [t,z,y,x] rank values'
351        quit(-1)
352
353    return cape, cin, zlfc, plfc, li, psldims, pslvdims
354
355
356
357def var_clt(cfra):
358    """ Function to compute the total cloud fraction following 'newmicro.F90' from
359      LMDZ using 1D vertical column values
360      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
361    """
362    ZEPSEC=1.0E-12
363
364    fname = 'var_clt'
365
366    zclear = 1.
367    zcloud = 0.
368
369    dz = cfra.shape[0]
370    for iz in range(dz):
371        zclear =zclear*(1.-np.max([cfra[iz],zcloud]))/(1.-np.min([zcloud,1.-ZEPSEC]))
372        clt = 1. - zclear
373        zcloud = cfra[iz]
374
375    return clt
376
377def compute_clt(cldfra, dimns, dimvns):
378    """ Function to compute the total cloud fraction following 'newmicro.F90' from
379      LMDZ
380    compute_clt(cldfra, dimnames)
381      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
382      [dimns]= list of the name of the dimensions of [cldfra]
383      [dimvns]= list of the name of the variables with the values of the
384        dimensions of [cldfra]
385    """
386    fname = 'compute_clt'
387
388    cltdims = dimns[:]
389    cltvdims = dimvns[:]
390
391    if len(cldfra.shape) == 4:
392        clt = np.zeros((cldfra.shape[0],cldfra.shape[2],cldfra.shape[3]),            \
393          dtype=np.float)
394        dx = cldfra.shape[3]
395        dy = cldfra.shape[2]
396        dz = cldfra.shape[1]
397        dt = cldfra.shape[0]
398        cltdims.pop(1)
399        cltvdims.pop(1)
400
401        for it in range(dt):
402            for ix in range(dx):
403                for iy in range(dy):
404                    zclear = 1.
405                    zcloud = 0.
406                    gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
407                    clt[it,iy,ix] = var_clt(cldfra[it,:,iy,ix])
408
409    else:
410        clt = np.zeros((cldfra.shape[1],cldfra.shape[2]), dtype=np.float)
411        dx = cldfra.shape[2]
412        dy = cldfra.shape[1]
413        dy = cldfra.shape[0]
414        cltdims.pop(0)
415        cltvdims.pop(0)
416        for ix in range(dx):
417            for iy in range(dy):
418                zclear = 1.
419                zcloud = 0.
420                gen.percendone(ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
421                clt[iy,ix] = var_clt(cldfra[:,iy,ix])
422
423    return clt, cltdims, cltvdims
424
425def Forcompute_clt(cldfra, dimns, dimvns):
426    """ Function to compute the total cloud fraction following 'newmicro.F90' from
427      LMDZ via a Fortran module
428    compute_clt(cldfra, dimnames)
429      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
430      [dimns]= list of the name of the dimensions of [cldfra]
431      [dimvns]= list of the name of the variables with the values of the
432        dimensions of [cldfra]
433    """
434    fname = 'Forcompute_clt'
435
436    cltdims = dimns[:]
437    cltvdims = dimvns[:]
438
439
440    if len(cldfra.shape) == 4:
441        clt = np.zeros((cldfra.shape[0],cldfra.shape[2],cldfra.shape[3]),            \
442          dtype=np.float)
443        dx = cldfra.shape[3]
444        dy = cldfra.shape[2]
445        dz = cldfra.shape[1]
446        dt = cldfra.shape[0]
447        cltdims.pop(1)
448        cltvdims.pop(1)
449
450        clt = fdin.module_fordiagnostics.compute_clt4d2(cldfra[:])
451
452    else:
453        clt = np.zeros((cldfra.shape[1],cldfra.shape[2]), dtype=np.float)
454        dx = cldfra.shape[2]
455        dy = cldfra.shape[1]
456        dy = cldfra.shape[0]
457        cltdims.pop(0)
458        cltvdims.pop(0)
459
460        clt = fdin.module_fordiagnostics.compute_clt3d1(cldfra[:])
461
462    return clt, cltdims, cltvdims
463
464def var_cllmh(cfra, p):
465    """ Fcuntion to compute cllmh on a 1D column
466    """
467
468    fname = 'var_cllmh'
469
470    ZEPSEC =1.0E-12
471    prmhc = 440.*100.
472    prmlc = 680.*100.
473
474    zclearl = 1.
475    zcloudl = 0.
476    zclearm = 1.
477    zcloudm = 0.
478    zclearh = 1.
479    zcloudh = 0.
480
481    dvz = cfra.shape[0]
482
483    cllmh = np.ones((3), dtype=np.float)
484
485    for iz in range(dvz):
486        if p[iz] < prmhc:
487            cllmh[2] = cllmh[2]*(1.-np.max([cfra[iz], zcloudh]))/(1.-                \
488              np.min([zcloudh,1.-ZEPSEC]))
489            zcloudh = cfra[iz]
490        elif p[iz] >= prmhc and p[iz] < prmlc:
491            cllmh[1] = cllmh[1]*(1.-np.max([cfra[iz], zcloudm]))/(1.-                \
492              np.min([zcloudm,1.-ZEPSEC]))
493            zcloudm = cfra[iz]
494        elif p[iz] >= prmlc:
495            cllmh[0] = cllmh[0]*(1.-np.max([cfra[iz], zcloudl]))/(1.-                \
496              np.min([zcloudl,1.-ZEPSEC]))
497            zcloudl = cfra[iz]
498
499    cllmh = 1.- cllmh
500
501    return cllmh
502
503def Forcompute_cllmh(cldfra, pres, dimns, dimvns):
504    """ Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ via Fortran subroutine
505    compute_clt(cldfra, pres, dimns, dimvns)
506      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
507      [pres] = pressure field
508      [dimns]= list of the name of the dimensions of [cldfra]
509      [dimvns]= list of the name of the variables with the values of the
510        dimensions of [cldfra]
511    """
512    fname = 'Forcompute_cllmh'
513
514    cllmhdims = dimns[:]
515    cllmhvdims = dimvns[:]
516
517    if len(cldfra.shape) == 4:
518        dx = cldfra.shape[3]
519        dy = cldfra.shape[2]
520        dz = cldfra.shape[1]
521        dt = cldfra.shape[0]
522        cllmhdims.pop(1)
523        cllmhvdims.pop(1)
524
525        cllmh = fdin.module_fordiagnostics.compute_cllmh4d2(cldfra[:], pres[:])
526       
527    else:
528        dx = cldfra.shape[2]
529        dy = cldfra.shape[1]
530        dz = cldfra.shape[0]
531        cllmhdims.pop(0)
532        cllmhvdims.pop(0)
533
534        cllmh = fdin.module_fordiagnostics.compute_cllmh3d1(cldfra[:], pres[:])
535
536    return cllmh, cllmhdims, cllmhvdims
537
538def compute_cllmh(cldfra, pres, dimns, dimvns):
539    """ Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ
540    compute_clt(cldfra, pres, dimns, dimvns)
541      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
542      [pres] = pressure field
543      [dimns]= list of the name of the dimensions of [cldfra]
544      [dimvns]= list of the name of the variables with the values of the
545        dimensions of [cldfra]
546    """
547    fname = 'compute_cllmh'
548
549    cllmhdims = dimns[:]
550    cllmhvdims = dimvns[:]
551
552    if len(cldfra.shape) == 4:
553        dx = cldfra.shape[3]
554        dy = cldfra.shape[2]
555        dz = cldfra.shape[1]
556        dt = cldfra.shape[0]
557        cllmhdims.pop(1)
558        cllmhvdims.pop(1)
559
560        cllmh = np.ones(tuple([3, dt, dy, dx]), dtype=np.float)
561
562        for it in range(dt):
563            for ix in range(dx):
564                for iy in range(dy):
565                    gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
566                    cllmh[:,it,iy,ix] = var_cllmh(cldfra[it,:,iy,ix], pres[it,:,iy,ix])
567       
568    else:
569        dx = cldfra.shape[2]
570        dy = cldfra.shape[1]
571        dz = cldfra.shape[0]
572        cllmhdims.pop(0)
573        cllmhvdims.pop(0)
574
575        cllmh = np.ones(tuple([3, dy, dx]), dtype=np.float)
576
577        for ix in range(dx):
578            for iy in range(dy):
579                gen.percendone(ix*dy + iy,dx*dy, 5, 'diagnosted')
580                cllmh[:,iy,ix] = var_cllmh(cldfra[:,iy,ix], pres[:,iy,ix])
581
582    return cllmh, cllmhdims, cllmhvdims
583
584def compute_clivi(dens, qtot, dimns, dimvns):
585    """ Function to compute cloud-ice water path (clivi)
586      [dens] = density [in kgkg-1] (assuming [t],z,y,x)
587      [qtot] = added mixing ratio of all cloud-ice species in [kgkg-1] (assuming [t],z,y,x)
588      [dimns]= list of the name of the dimensions of [q]
589      [dimvns]= list of the name of the variables with the values of the
590        dimensions of [q]
591    """
592    fname = 'compute_clivi'
593
594    clividims = dimns[:]
595    clivivdims = dimvns[:]
596
597    if len(qtot.shape) == 4:
598        clividims.pop(1)
599        clivivdims.pop(1)
600    else:
601        clividims.pop(0)
602        clivivdims.pop(0)
603
604    data1 = dens*qtot
605    clivi = np.sum(data1, axis=1)
606
607    return clivi, clividims, clivivdims
608
609
610def compute_clwvl(dens, qtot, dimns, dimvns):
611    """ Function to compute condensed water path (clwvl)
612      [dens] = density [in kgkg-1] (assuming [t],z,y,x)
613      [qtot] = added mixing ratio of all cloud-water species in [kgkg-1] (assuming [t],z,y,x)
614      [dimns]= list of the name of the dimensions of [q]
615      [dimvns]= list of the name of the variables with the values of the
616        dimensions of [q]
617    """
618    fname = 'compute_clwvl'
619
620    clwvldims = dimns[:]
621    clwvlvdims = dimvns[:]
622
623    if len(qtot.shape) == 4:
624        clwvldims.pop(1)
625        clwvlvdims.pop(1)
626    else:
627        clwvldims.pop(0)
628        clwvlvdims.pop(0)
629
630    data1 = dens*qtot
631    clwvl = np.sum(data1, axis=1)
632
633    return clwvl, clwvldims, clwvlvdims
634
635def var_virtualTemp (temp,rmix):
636    """ This function returns virtual temperature in K,
637      temp: temperature [K]
638      rmix: mixing ratio in [kgkg-1]
639    """
640
641    fname = 'var_virtualTemp'
642
643    virtual=temp*(0.622+rmix)/(0.622*(1.+rmix))
644
645    return virtual
646
647def Forcompute_zint(var, zinterlev, zweights, dimns, dimvns):
648    """ Function to compute a vertical integration of volumetric quantities
649    Forcompute_mrso(smois, dsoil, dimns, dimvns)
650      [var]= values (assuming [[t],z,y,x]) [volumetric units]
651      [zinterlev]= depth of each layer (assuming [z]) [same z units as var]
652      [zweights]= weights to apply to each level (just in case...)
653      [dimns]= list of the name of the dimensions of [smois]
654      [dimvns]= list of the name of the variables with the values of the
655        dimensions of [smois]
656    """
657    fname = 'Forcompute_zint'
658
659    vardims = dimns[:]
660    varvdims = dimvns[:]
661
662    if len(var.shape) == 4:
663        zint = np.zeros((var.shape[0],var.shape[2],var.shape[3]), dtype=np.float)
664        dx = var.shape[3]
665        dy = var.shape[2]
666        dz = var.shape[1]
667        dt = var.shape[0]
668        vardims.pop(1)
669        varvdims.pop(1)
670
671        zintvart=fdin.module_fordiagnostics.compute_zint4d(var4d=var[:].transpose(), \
672          dlev=zinterlev[:].transpose(), zweight=zweights[:].transpose(), d1=dx,     \
673          d2=dy, d3=dz, d4=dt)
674        zintvar = zintvart.transpose()
675    else:
676        print errormsg
677        print '  ' + fname + ': rank', len(var.shape), 'not ready !!'
678        print '  it only computes 4D [t,z,y,x] rank values'
679        quit(-1)
680
681    return zintvar, vardims, varvdims
682
683def var_mslp(pres, psfc, ter, tk, qv):
684    """ Function to compute mslp on a 1D column
685    """
686
687    fname = 'var_mslp'
688
689    N = 1.0
690    expon=287.04*.0065/9.81
691    pref = 40000.
692
693# First find where about 400 hPa is located
694    dz=len(pres) 
695
696    kref = -1
697    pinc = pres[0] - pres[dz-1]
698
699    if pinc < 0.:
700        for iz in range(1,dz):
701            if pres[iz-1] >= pref and pres[iz] < pref: 
702                kref = iz
703                break
704    else:
705        for iz in range(dz-1):
706            if pres[iz] >= pref and pres[iz+1] < pref: 
707                kref = iz
708                break
709
710    if kref == -1:
711        print errormsg
712        print '  ' + fname + ': no reference pressure:',pref,'found!!'
713        print '    values:',pres[:]
714        quit(-1)
715
716    mslp = 0.
717
718# We are below both the ground and the lowest data level.
719
720# First, find the model level that is closest to a "target" pressure
721# level, where the "target" pressure is delta-p less that the local
722# value of a horizontally smoothed surface pressure field.  We use
723# delta-p = 150 hPa here. A standard lapse rate temperature profile
724# passing through the temperature at this model level will be used
725# to define the temperature profile below ground.  This is similar
726# to the Benjamin and Miller (1990) method, using 
727# 700 hPa everywhere for the "target" pressure.
728
729# ptarget = psfc - 15000.
730    ptarget = 70000.
731    dpmin=1.e4
732    kupper = 0
733    if pinc > 0.:
734        for iz in range(dz-1,0,-1):
735            kupper = iz
736            dp=np.abs( pres[iz] - ptarget )
737            if dp < dpmin: exit
738            dpmin = np.min([dpmin, dp])
739    else:
740        for iz in range(dz):
741            kupper = iz
742            dp=np.abs( pres[iz] - ptarget )
743            if dp < dpmin: exit
744            dpmin = np.min([dpmin, dp])
745
746    pbot=np.max([pres[0], psfc])
747#    zbot=0.
748
749#    tbotextrap=tk(i,j,kupper,itt)*(pbot/pres_field(i,j,kupper,itt))**expon
750#    tvbotextrap=virtual(tbotextrap,qv(i,j,1,itt))
751
752#    data_out(i,j,itt,1) = (zbot+tvbotextrap/.0065*(1.-(interp_levels(1)/pbot)**expon))
753    tbotextrap = tk[kupper]*(psfc/ptarget)**expon
754    tvbotextrap = var_virtualTemp(tbotextrap, qv[kupper])
755    mslp = psfc*( (tvbotextrap+0.0065*ter)/tvbotextrap)**(1./expon)
756
757    return mslp
758
759def compute_mslp(pressure, psurface, terrain, temperature, qvapor, dimns, dimvns):
760    """ Function to compute mslp: mean sea level pressure following p_interp.F90 from WRF
761    var_mslp(pres, ter, tk, qv, dimns, dimvns)
762      [pressure]= pressure field [Pa] (assuming [[t],z,y,x])
763      [psurface]= surface pressure field [Pa]
764      [terrain]= topography [m]
765      [temperature]= temperature [K]
766      [qvapor]= water vapour mixing ratio [kgkg-1]
767      [dimns]= list of the name of the dimensions of [cldfra]
768      [dimvns]= list of the name of the variables with the values of the
769        dimensions of [pres]
770    """
771
772    fname = 'compute_mslp'
773
774    mslpdims = list(dimns[:])
775    mslpvdims = list(dimvns[:])
776
777    if len(pressure.shape) == 4:
778        mslpdims.pop(1)
779        mslpvdims.pop(1)
780    else:
781        mslpdims.pop(0)
782        mslpvdims.pop(0)
783
784    if len(pressure.shape) == 4:
785        dx = pressure.shape[3]
786        dy = pressure.shape[2]
787        dz = pressure.shape[1]
788        dt = pressure.shape[0]
789
790        mslpv = np.zeros(tuple([dt, dy, dx]), dtype=np.float)
791
792# Terrain... to 2D !
793        terval = np.zeros(tuple([dy, dx]), dtype=np.float)
794        if len(terrain.shape) == 3:
795            terval = terrain[0,:,:]
796        else:
797            terval = terrain
798
799        for ix in range(dx):
800            for iy in range(dy):
801                if terval[iy,ix] > 0.:
802                    for it in range(dt):
803                        mslpv[it,iy,ix] = var_mslp(pressure[it,:,iy,ix],             \
804                          psurface[it,iy,ix], terval[iy,ix], temperature[it,:,iy,ix],\
805                          qvapor[it,:,iy,ix])
806
807                        gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
808                else:
809                    mslpv[:,iy,ix] = psurface[:,iy,ix]
810
811    else:
812        dx = pressure.shape[2]
813        dy = pressure.shape[1]
814        dz = pressure.shape[0]
815
816        mslpv = np.zeros(tuple([dy, dx]), dtype=np.float)
817
818# Terrain... to 2D !
819        terval = np.zeros(tuple([dy, dx]), dtype=np.float)
820        if len(terrain.shape) == 3:
821            terval = terrain[0,:,:]
822        else:
823            terval = terrain
824
825        for ix in range(dx):
826            for iy in range(dy):
827                gen.percendone(ix*dy + iy,dx*dy, 5, 'diagnosted')
828                if terval[iy,ix] > 0.:
829                    mslpv[iy,ix] = var_mslp(pressure[:,iy,ix], psurface[iy,ix],      \
830                      terval[iy,ix], temperature[:,iy,ix], qvapor[:,iy,ix])
831                else:
832                    mslpv[iy,ix] = psfc[iy,ix]
833
834    return mslpv, mslpdims, mslpvdims
835
836def Forcompute_psl_ecmwf(ps, hgt, ta1, pa2, unpa1, dimns, dimvns):
837    """ Function to compute the sea-level pressure following Mats Hamrud and Philippe Courtier [Pa]
838    Forcompute_psl_ptarget(ps, hgt, ta1, pa2, unpa1, dimns, dimvns)
839      [ps]= surface pressure values (assuming [[t],y,x]) [Pa]
840      [hgt]= opography (assuming [y,x]) [m]
841      [ta1]= air-temperature values at first half-mass level (assuming [[t],y,x]) [K]
842      [pa2]= pressure values at second full-mass levels (assuming [[t],y,x]) [Pa]
843      [unpa1]= pressure values at first half-mass levels (assuming [[t],y,x]) [Pa]
844      [dimns]= list of the name of the dimensions of [pa]
845      [dimvns]= list of the name of the variables with the values of the
846        dimensions of [pa]
847    """
848    fname = 'Forcompute_psl_ecmwf'
849
850    vardims = dimns[:]
851    varvdims = dimvns[:]
852
853    if len(pa2.shape) == 3:
854        psl = np.zeros((pa2.shape[0],pa2.shape[1],pa2.shape[2]), dtype=np.float)
855        dx = pa2.shape[2]
856        dy = pa2.shape[1]
857        dt = pa2.shape[0]
858        pslt= fdin.module_fordiagnostics.compute_psl_ecmwf( ps=ps[:].transpose(),    \
859          hgt=hgt[:].transpose(), t=ta1[:].transpose(), press=pa2[:].transpose(),    \
860          unpress=unpa1[:].transpose(), d1=dx, d2=dy, d4=dt)
861        psl = pslt.transpose()
862    else:
863        print errormsg
864        print '  ' + fname + ': rank', len(pa2.shape), 'not ready !!'
865        print '  it only computes 3D [t,y,x] rank values'
866        quit(-1)
867
868    return psl, vardims, varvdims
869
870def Forcompute_psl_ptarget(pa, ps, ta, hgt, qv, target_pressure, dimns, dimvns):
871    """ Function to compute the sea-level pressure following target_pressure value
872      found in `p_interp.F'
873    Forcompute_psl_ptarget(pa, ps, ta, hgt, qv, dimns, dimvns)
874      [pa]= pressure values (assuming [[t],z,y,x]) [Pa]
875      [ps]= surface pressure values (assuming [[t],y,x]) [Pa]
876      [ta]= air-temperature values (assuming [[t],z,y,x]) [K]
877      [hgt]= opography (assuming [y,x]) [m]
878      [qv]= water vapour mixing ratio (assuming [[t],z,y,x]) [kgkg-1]
879      [dimns]= list of the name of the dimensions of [pa]
880      [dimvns]= list of the name of the variables with the values of the
881        dimensions of [pa]
882    """
883    fname = 'Forcompute_psl_ptarget'
884
885    psldims = dimns[:]
886    pslvdims = dimvns[:]
887
888    if len(pa.shape) == 4:
889        psl = np.zeros((pa.shape[0],pa.shape[2],pa.shape[3]), dtype=np.float)
890        dx = pa.shape[3]
891        dy = pa.shape[2]
892        dz = pa.shape[1]
893        dt = pa.shape[0]
894        psldims.pop(1)
895        pslvdims.pop(1)
896
897        pslt= fdin.module_fordiagnostics.compute_psl_ptarget4d2(                     \
898          press=pa[:].transpose(), ps=ps[:].transpose(), hgt=hgt[:].transpose(),     \
899          ta=ta[:].transpose(), qv=qv[:].transpose(), ptarget=target_pressure,       \
900          d1=dx, d2=dy, d3=dz, d4=dt)
901        psl = pslt.transpose()
902    else:
903        print errormsg
904        print '  ' + fname + ': rank', len(pa.shape), 'not ready !!'
905        print '  it only computes 4D [t,z,y,x] rank values'
906        quit(-1)
907
908    return psl, psldims, pslvdims
909
910def Forcompute_zmla_gen(theta, qratio, zpl, hgt, dimns, dimvns):
911    """ Function to compute the boundary layer height following a generic method
912      with Fortran
913    Forcompute_zmla_gen(theta, qratio, zpl, hgt, zmla, dimns, dimvns)
914      [theta]= potential air-temperature values (assuming [[t],z,y,x]) [K]
915      [qratio]= water mixing ratio (assuming [[t],z,y,x]) [kgkg-1]
916      [zpl]= height from sea level (assuming [[t],z,y,x]) [m]
917      [hgt]= topographical height (assuming [m]
918      [zmla]= boundary layer height [m]
919      [dimns]= list of the name of the dimensions of [theta]
920      [dimvns]= list of the name of the variables with the values of the
921        dimensions of [theta]
922    """
923    fname = 'Forcompute_zmla_gen'
924
925    zmladims = dimns[:]
926    zmlavdims = dimvns[:]
927
928    if len(theta.shape) == 4:
929        zmla= np.zeros((theta.shape[0],theta.shape[2],theta.shape[3]), dtype=np.float)
930
931        dx = theta.shape[3]
932        dy = theta.shape[2]
933        dz = theta.shape[1]
934        dt = theta.shape[0]
935        zmladims.pop(1)
936        zmlavdims.pop(1)
937
938        pzmla= fdin.module_fordiagnostics.compute_zmla_generic4d(                    \
939          tpot=theta[:].transpose(), qratio=qratio[:].transpose(),                   \
940          z=zpl[:].transpose(), hgt=hgt.transpose(), d1=dx, d2=dy, d3=dz, d4=dt)
941        zmla = pzmla.transpose()
942    else:
943        print errormsg
944        print '  ' + fname + ': rank', len(theta.shape), 'not ready !!'
945        print '  it only computes 4D [t,z,y,x] rank values'
946        quit(-1)
947
948    return zmla, zmladims, zmlavdims
949
950def Forcompute_zwind(ua, va, z, uas, vas, sina, cosa, zval, dimns, dimvns):
951    """ Function to compute the wind at a given height following the power law method
952    Forcompute_zwind(ua, va, zsl, uas, vas, hgt, sina, cosa, zval, dimns, dimvns)
953      [ua]= x-component of unstaggered 3D wind (assuming [[t],z,y,x]) [ms-1]
954      [va]= y-component of unstaggered 3D wind (assuming [[t],z,y,x]) [ms-1]
955      [z]= height above surface [m]
956      [uas]= x-component of unstaggered 10 m wind (assuming [[t],z,y,x]) [ms-1]
957      [vas]= y-component of unstaggered 10 m wind (assuming [[t],z,y,x]) [ms-1]
958      [sina]= local sine of map rotation [1.]
959      [cosa]= local cosine of map rotation [1.]
960      [zval]= desired height for winds [m]
961      [dimns]= list of the name of the dimensions of [ua]
962      [dimvns]= list of the name of the variables with the values of the
963        dimensions of [ua]
964    """
965    fname = 'Forcompute_zwind'
966
967    vardims = dimns[:]
968    varvdims = dimvns[:]
969
970    if len(ua.shape) == 4:
971        var1= np.zeros((ua.shape[0],ua.shape[2],ua.shape[3]), dtype=np.float)
972        var2= np.zeros((ua.shape[0],ua.shape[2],ua.shape[3]), dtype=np.float)
973
974        dx = ua.shape[3]
975        dy = ua.shape[2]
976        dz = ua.shape[1]
977        dt = ua.shape[0]
978        vardims.pop(1)
979        varvdims.pop(1)
980
981        pvar1, pvar2= fdin.module_fordiagnostics.compute_zwind4d(ua=ua.transpose(),  \
982          va=va[:].transpose(), z=z[:].transpose(), uas=uas.transpose(),             \
983          vas=vas.transpose(), sina=sina.transpose(), cosa=cosa.transpose(),         \
984          zextrap=zval, d1=dx, d2=dy, d3=dz, d4=dt)
985        var1 = pvar1.transpose()
986        var2 = pvar2.transpose()
987    else:
988        print errormsg
989        print '  ' + fname + ': rank', len(ua.shape), 'not ready !!'
990        print '  it only computes 4D [t,z,y,x] rank values'
991        quit(-1)
992
993    return var1, var2, vardims, varvdims
994
995def Forcompute_zwind_log(ua, va, z, uas, vas, sina, cosa, zval, dimns, dimvns):
996    """ Function to compute the wind at a given height following the logarithmic law method
997    Forcompute_zwind(ua, va, zsl, uas, vas, hgt, sina, cosa, zval, dimns, dimvns)
998      [ua]= x-component of unstaggered 3D wind (assuming [[t],z,y,x]) [ms-1]
999      [va]= y-component of unstaggered 3D wind (assuming [[t],z,y,x]) [ms-1]
1000      [z]= height above surface [m]
1001      [uas]= x-component of unstaggered 10 m wind (assuming [[t],z,y,x]) [ms-1]
1002      [vas]= y-component of unstaggered 10 m wind (assuming [[t],z,y,x]) [ms-1]
1003      [sina]= local sine of map rotation [1.]
1004      [cosa]= local cosine of map rotation [1.]
1005      [zval]= desired height for winds [m]
1006      [dimns]= list of the name of the dimensions of [ua]
1007      [dimvns]= list of the name of the variables with the values of the
1008        dimensions of [ua]
1009    """
1010    fname = 'Forcompute_zwind_log'
1011
1012    vardims = dimns[:]
1013    varvdims = dimvns[:]
1014
1015    if len(ua.shape) == 4:
1016        var1= np.zeros((ua.shape[0],ua.shape[2],ua.shape[3]), dtype=np.float)
1017        var2= np.zeros((ua.shape[0],ua.shape[2],ua.shape[3]), dtype=np.float)
1018
1019        dx = ua.shape[3]
1020        dy = ua.shape[2]
1021        dz = ua.shape[1]
1022        dt = ua.shape[0]
1023        vardims.pop(1)
1024        varvdims.pop(1)
1025
1026        pvar1, pvar2= fdin.module_fordiagnostics.compute_zwind_log4d(                \
1027          ua=ua.transpose(), va=va[:].transpose(), z=z[:].transpose(),               \
1028          uas=uas.transpose(), vas=vas.transpose(), sina=sina.transpose(),           \
1029          cosa=cosa.transpose(), zextrap=zval, d1=dx, d2=dy, d3=dz, d4=dt)
1030        var1 = pvar1.transpose()
1031        var2 = pvar2.transpose()
1032    else:
1033        print errormsg
1034        print '  ' + fname + ': rank', len(ua.shape), 'not ready !!'
1035        print '  it only computes 4D [t,z,y,x] rank values'
1036        quit(-1)
1037
1038    return var1, var2, vardims, varvdims
1039
1040def Forcompute_zwindMO(ust, znt, rmol, uas, vas, sina, cosa, zval, dimns, dimvns):
1041    """ Function to compute the wind at a given height following the Monin-Obukhov theory
1042    Forcompute_zwind(ust, znt, rmol, uas, vas, sina, cosa, zval, dimns, dimvns)
1043      [ust]: u* in similarity theory (assuming [[t],y,x]) [ms-1]
1044      [znt]: thermal time-varying roughness length (assuming [[t],y,x]) [m]
1045      [rmol]: inverse of Obukhov length (assuming [[t],y,x]) [m-1]
1046      [uas]= x-component of unstaggered 10 m wind (assuming [[t],y,x]) [ms-1]
1047      [vas]= y-component of unstaggered 10 m wind (assuming [[t],y,x]) [ms-1]
1048      [sina]= local sine of map rotation [1.]
1049      [cosa]= local cosine of map rotation [1.]
1050      [zval]= desired height for winds [m]
1051      [dimns]= list of the name of the dimensions of [uas]
1052      [dimvns]= list of the name of the variables with the values of the
1053        dimensions of [uas]
1054    """
1055    fname = 'Forcompute_zwindMO'
1056
1057    vardims = dimns[:]
1058    varvdims = dimvns[:]
1059
1060    if len(uas.shape) == 3:
1061        var1= np.zeros((uas.shape[0],uas.shape[1],uas.shape[2]), dtype=np.float)
1062        var2= np.zeros((uas.shape[0],uas.shape[1],uas.shape[2]), dtype=np.float)
1063
1064        dx = uas.shape[2]
1065        dy = uas.shape[1]
1066        dt = uas.shape[0]
1067
1068        pvar1, pvar2 = fdin.module_fordiagnostics.compute_zwindmo3d(                 \
1069          ust=ust.transpose(), znt=znt[:].transpose(), rmol=rmol[:].transpose(),     \
1070          uas=uas.transpose(), vas=vas.transpose(), sina=sina.transpose(),           \
1071          cosa=cosa.transpose(), newz=zval, d1=dx, d2=dy, d3=dt)
1072        var1 = pvar1.transpose()
1073        var2 = pvar2.transpose()
1074    else:
1075        print errormsg
1076        print '  ' + fname + ': rank', len(uas.shape), 'not ready !!'
1077        print '  it only computes 3D [t,y,x] rank values'
1078        quit(-1)
1079
1080    return var1, var2, vardims, varvdims
1081
1082def compute_OMEGAw(omega, p, t, dimns, dimvns):
1083    """ Function to transform OMEGA [Pas-1] to velocities [ms-1]
1084      tacking: https://www.ncl.ucar.edu/Document/Functions/Contributed/omega_to_w.shtml
1085      [omega] = vertical velocity [in ms-1] (assuming [t],z,y,x)
1086      [p] = pressure in [Pa] (assuming [t],z,y,x)
1087      [t] = temperature in [K] (assuming [t],z,y,x)
1088      [dimns]= list of the name of the dimensions of [q]
1089      [dimvns]= list of the name of the variables with the values of the
1090        dimensions of [q]
1091    """
1092    fname = 'compute_OMEGAw'
1093
1094    rgas = 287.058            # J/(kg-K) => m2/(s2 K)
1095    g    = 9.80665            # m/s2
1096
1097    wdims = dimns[:]
1098    wvdims = dimvns[:]
1099
1100    rho  = p/(rgas*t)         # density => kg/m3
1101    w    = -omega/(rho*g)     
1102
1103    return w, wdims, wvdims
1104
1105def compute_prw(dens, q, dimns, dimvns):
1106    """ Function to compute water vapour path (prw)
1107      [dens] = density [in kgkg-1] (assuming [t],z,y,x)
1108      [q] = mixing ratio in [kgkg-1] (assuming [t],z,y,x)
1109      [dimns]= list of the name of the dimensions of [q]
1110      [dimvns]= list of the name of the variables with the values of the
1111        dimensions of [q]
1112    """
1113    fname = 'compute_prw'
1114
1115    prwdims = dimns[:]
1116    prwvdims = dimvns[:]
1117
1118    if len(q.shape) == 4:
1119        prwdims.pop(1)
1120        prwvdims.pop(1)
1121    else:
1122        prwdims.pop(0)
1123        prwvdims.pop(0)
1124
1125    data1 = dens*q
1126    prw = np.sum(data1, axis=1)
1127
1128    return prw, prwdims, prwvdims
1129
1130def compute_rh(p, t, q, dimns, dimvns):
1131    """ Function to compute relative humidity following 'Tetens' equation (T,P) ...'
1132      [t]= temperature (assuming [[t],z,y,x] in [K])
1133      [p] = pressure field (assuming in [hPa])
1134      [q] = mixing ratio in [kgkg-1]
1135      [dimns]= list of the name of the dimensions of [t]
1136      [dimvns]= list of the name of the variables with the values of the
1137        dimensions of [t]
1138    """
1139    fname = 'compute_rh'
1140
1141    rhdims = dimns[:]
1142    rhvdims = dimvns[:]
1143
1144    data1 = 10.*0.6112*np.exp(17.67*(t-273.16)/(t-29.65))
1145    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
1146
1147    rh = q/data2
1148
1149    return rh, rhdims, rhvdims
1150
1151def compute_td(p, temp, qv, dimns, dimvns):
1152    """ Function to compute the dew point temperature
1153      [p]= pressure [Pa]
1154      [temp]= temperature [C]
1155      [qv]= mixing ratio [kgkg-1]
1156      [dimns]= list of the name of the dimensions of [p]
1157      [dimvns]= list of the name of the variables with the values of the
1158        dimensions of [p]
1159    """
1160    fname = 'compute_td'
1161
1162#    print '    ' + fname + ': computing dew-point temperature from TS as t and Tetens...'
1163# tacking from: http://en.wikipedia.org/wiki/Dew_point
1164    tk = temp
1165    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
1166    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
1167
1168    rh = qv/data2
1169               
1170    pa = rh * data1
1171    td = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
1172
1173    tddims = dimns[:]
1174    tdvdims = dimvns[:]
1175
1176    return td, tddims, tdvdims
1177
1178def var_WRFtime(timewrfv, refdate='19491201000000', tunitsval='minutes'):
1179    """ Function to copmute CFtimes from WRFtime variable
1180      refdate= [YYYYMMDDMIHHSS] format of reference date
1181      tunitsval= CF time units
1182      timewrfv= matrix string values of WRF 'Times' variable
1183    """
1184    fname = 'var_WRFtime'
1185
1186    yrref=refdate[0:4]
1187    monref=refdate[4:6]
1188    dayref=refdate[6:8]
1189    horref=refdate[8:10]
1190    minref=refdate[10:12]
1191    secref=refdate[12:14]
1192
1193    refdateS = yrref + '-' + monref + '-' + dayref + ' ' + horref + ':' + minref +   \
1194      ':' + secref
1195
1196    dt = timewrfv.shape[0]
1197    WRFtime = np.zeros((dt), dtype=np.float)
1198
1199    for it in range(dt):
1200        wrfdates = gen.datetimeStr_conversion(timewrfv[it,:],'WRFdatetime', 'matYmdHMS')
1201        WRFtime[it] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
1202
1203    tunits = tunitsval + ' since ' + refdateS
1204
1205    return WRFtime, tunits
1206
1207def turbulence_var(varv, dimvn, dimn):
1208    """ Function to compute the Taylor's decomposition turbulence term from a a given variable
1209      x*=<x^2>_t-(<X>_t)^2
1210    turbulence_var(varv,dimn)
1211      varv= values of the variable
1212      dimvn= names of the dimension of the variable
1213      dimn= names of the dimensions (as a dictionary with 'X', 'Y', 'Z', 'T')
1214    >>> turbulence_var(np.arange((27)).reshape(3,3,3),['time','y','x'],{'T':'time', 'Y':'y', 'X':'x'})
1215    [[ 54.  54.  54.]
1216     [ 54.  54.  54.]
1217     [ 54.  54.  54.]]
1218    """
1219    fname = 'turbulence_varv'
1220
1221    timedimid = dimvn.index(dimn['T'])
1222
1223    varv2 = varv*varv
1224
1225    vartmean = np.mean(varv, axis=timedimid)
1226    var2tmean = np.mean(varv2, axis=timedimid)
1227
1228    varvturb = var2tmean - (vartmean*vartmean)
1229
1230    return varvturb
1231
1232def compute_turbulence(v, dimns, dimvns):
1233    """ Function to compute the rubulence term of the Taylor's decomposition ...'
1234      x*=<x^2>_t-(<X>_t)^2
1235      [v]= variable (assuming [[t],z,y,x])
1236      [dimns]= list of the name of the dimensions of [v]
1237      [dimvns]= list of the name of the variables with the values of the
1238        dimensions of [v]
1239    """
1240    fname = 'compute_turbulence'
1241
1242    turbdims = dimns[:]
1243    turbvdims = dimvns[:]
1244
1245    turbdims.pop(0)
1246    turbvdims.pop(0)
1247
1248    v2 = v*v
1249
1250    vartmean = np.mean(v, axis=0)
1251    var2tmean = np.mean(v2, axis=0)
1252
1253    turb = var2tmean - (vartmean*vartmean)
1254
1255    return turb, turbdims, turbvdims
1256
1257def compute_wds(u, v, dimns, dimvns):
1258    """ Function to compute the wind direction
1259      [u]= W-E wind direction [ms-1, knot, ...]
1260      [v]= N-S wind direction [ms-1, knot, ...]
1261      [dimns]= list of the name of the dimensions of [u]
1262      [dimvns]= list of the name of the variables with the values of the
1263        dimensions of [u]
1264    """
1265    fname = 'compute_wds'
1266
1267#    print '    ' + fname + ': computing wind direction as ATAN2(v,u) ...'
1268    theta = np.arctan2(v,u)
1269    theta = np.where(theta < 0., theta + 2.*np.pi, theta)
1270
1271    wds = 360.*theta/(2.*np.pi)
1272
1273    wdsdims = dimns[:]
1274    wdsvdims = dimvns[:]
1275
1276    return wds, wdsdims, wdsvdims
1277
1278def compute_wss(u, v, dimns, dimvns):
1279    """ Function to compute the wind speed
1280      [u]= W-E wind direction [ms-1, knot, ...]
1281      [v]= N-S wind direction [ms-1, knot, ...]
1282      [dimns]= list of the name of the dimensions of [u]
1283      [dimvns]= list of the name of the variables with the values of the
1284        dimensions of [u]
1285    """
1286    fname = 'compute_wss'
1287
1288#    print '    ' + fname + ': computing wind speed as SQRT(v**2 + u**2) ...'
1289    wss = np.sqrt(u*u + v*v)
1290
1291    wssdims = dimns[:]
1292    wssvdims = dimvns[:]
1293
1294    return wss, wssdims, wssvdims
1295
1296def timeunits_seconds(dtu):
1297    """ Function to transform a time units to seconds
1298    timeunits_seconds(timeuv)
1299      [dtu]= time units value to transform in seconds
1300    """
1301    fname='timunits_seconds'
1302
1303    if dtu == 'years':
1304        times = 365.*24.*3600.
1305    elif dtu == 'weeks':
1306        times = 7.*24.*3600.
1307    elif dtu == 'days':
1308        times = 24.*3600.
1309    elif dtu == 'hours':
1310        times = 3600.
1311    elif dtu == 'minutes':
1312        times = 60.
1313    elif dtu == 'seconds':
1314        times = 1.
1315    elif dtu == 'miliseconds':
1316        times = 1./1000.
1317    else:
1318        print errormsg
1319        print '  ' + fname  + ": time units '" + dtu + "' not ready !!"
1320        quit(-1)
1321
1322    return times
1323
1324def compute_WRFhur(t, p, qv, dimns, dimvns):
1325    """ Function to compute WRF relative humidity following Teten's equation
1326      t= orginal WRF temperature
1327      p= original WRF pressure (P + PB)
1328      formula:
1329          temp = theta*(p/p0)**(R/Cp)
1330
1331    """
1332    fname = 'compute_WRFtd'
1333
1334    tk = (t+300.)*(p/fdef.module_definitions.p0ref)**(fdef.module_definitions.rcp)
1335    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
1336    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
1337
1338    rh = qv/data2
1339               
1340    dnamesvar = ['Time','bottom_top','south_north','west_east']
1341    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1342
1343    return rh, dnamesvar, dvnamesvar
1344
1345def compute_WRFua(u, v, sina, cosa, dimns, dimvns):
1346    """ Function to compute geographical rotated WRF 3D winds
1347      u= orginal WRF x-wind
1348      v= orginal WRF y-wind
1349      sina= original WRF local sinus of map rotation
1350      cosa= original WRF local cosinus of map rotation
1351      formula:
1352        ua = u*cosa-va*sina
1353        va = u*sina+va*cosa
1354    """
1355    fname = 'compute_WRFua'
1356
1357    var0 = u
1358    var1 = v
1359    var2 = sina
1360    var3 = cosa
1361
1362    # un-staggering variables
1363    unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1364    ua = np.zeros(tuple(unstgdims), dtype=np.float)
1365    unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1366    unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1367    unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1368    unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1369
1370    for iz in range(var0.shape[1]):
1371        ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
1372
1373    dnamesvar = ['Time','bottom_top','south_north','west_east']
1374    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1375
1376    return ua, dnamesvar, dvnamesvar
1377
1378def compute_WRFva(u, v, sina, cosa, dimns, dimvns):
1379    """ Function to compute geographical rotated WRF 3D winds
1380      u= orginal WRF x-wind
1381      v= orginal WRF y-wind
1382      sina= original WRF local sinus of map rotation
1383      cosa= original WRF local cosinus of map rotation
1384      formula:
1385        ua = u*cosa-va*sina
1386        va = u*sina+va*cosa
1387    """
1388    fname = 'compute_WRFva'
1389
1390    var0 = u
1391    var1 = v
1392    var2 = sina
1393    var3 = cosa
1394
1395    # un-staggering variables
1396    unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1397    va = np.zeros(tuple(unstgdims), dtype=np.float)
1398    unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1399    unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1400    unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1401    unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1402
1403    for iz in range(var0.shape[1]):
1404        va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
1405
1406    dnamesvar = ['Time','bottom_top','south_north','west_east']
1407    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1408
1409    return va, dnamesvar, dvnamesvar
1410
1411def compute_WRFuava(u, v, sina, cosa, dimns, dimvns):
1412    """ Function to compute geographical rotated WRF 3D winds
1413      u= orginal WRF x-wind
1414      v= orginal WRF y-wind
1415      sina= original WRF local sinus of map rotation
1416      cosa= original WRF local cosinus of map rotation
1417      formula:
1418        ua = u*cosa-va*sina
1419        va = u*sina+va*cosa
1420    """
1421    fname = 'compute_WRFuava'
1422
1423    var0 = u
1424    var1 = v
1425    var2 = sina
1426    var3 = cosa
1427
1428    # un-staggering variables
1429    unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1430    ua = np.zeros(tuple(unstgdims), dtype=np.float)
1431    va = np.zeros(tuple(unstgdims), dtype=np.float)
1432    unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1433    unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1434    unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1435    unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1436
1437    for iz in range(var0.shape[1]):
1438        ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
1439        va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
1440
1441    dnamesvar = ['Time','bottom_top','south_north','west_east']
1442    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1443
1444    return ua, va, dnamesvar, dvnamesvar
1445
1446def compute_WRFuas(u10, v10, sina, cosa, dimns, dimvns):
1447    """ Function to compute geographical rotated WRF 2-meter x-wind
1448      u10= orginal WRF 10m x-wind
1449      v10= orginal WRF 10m y-wind
1450      sina= original WRF local sinus of map rotation
1451      cosa= original WRF local cosinus of map rotation
1452      formula:
1453        uas = u10*cosa-va10*sina
1454        vas = u10*sina+va10*cosa
1455    """
1456    fname = 'compute_WRFuas'
1457
1458    var0 = u10
1459    var1 = v10
1460    var2 = sina
1461    var3 = cosa
1462
1463    uas = np.zeros(var0.shape, dtype=np.float)
1464    vas = np.zeros(var0.shape, dtype=np.float)
1465
1466    uas = var0*var3 - var1*var2
1467
1468    dnamesvar = ['Time','south_north','west_east']
1469    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1470
1471    return uas, dnamesvar, dvnamesvar
1472
1473def compute_WRFvas(u10, v10, sina, cosa, dimns, dimvns):
1474    """ Function to compute geographical rotated WRF 2-meter y-wind
1475      u10= orginal WRF 10m x-wind
1476      v10= orginal WRF 10m y-wind
1477      sina= original WRF local sinus of map rotation
1478      cosa= original WRF local cosinus of map rotation
1479      formula:
1480        uas = u10*cosa-va10*sina
1481        vas = u10*sina+va10*cosa
1482    """
1483    fname = 'compute_WRFvas'
1484
1485    var0 = u10
1486    var1 = v10
1487    var2 = sina
1488    var3 = cosa
1489
1490    uas = np.zeros(var0.shape, dtype=np.float)
1491    vas = np.zeros(var0.shape, dtype=np.float)
1492
1493    vas = var0*var2 + var1*var3
1494
1495    dnamesvar = ['Time','south_north','west_east']
1496    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1497
1498    return vas, dnamesvar, dvnamesvar
1499
1500def compute_WRFuasvas(u10, v10, sina, cosa, dimns, dimvns):
1501    """ Function to compute geographical rotated WRF 2-meter winds
1502      u10= orginal WRF 10m x-wind
1503      v10= orginal WRF 10m y-wind
1504      sina= original WRF local sinus of map rotation
1505      cosa= original WRF local cosinus of map rotation
1506      formula:
1507        uas = u10*cosa-va10*sina
1508        vas = u10*sina+va10*cosa
1509    """
1510    fname = 'compute_WRFuasvas'
1511
1512    var0 = u10
1513    var1 = v10
1514    var2 = sina
1515    var3 = cosa
1516
1517    uas = np.zeros(var0.shape, dtype=np.float)
1518    vas = np.zeros(var0.shape, dtype=np.float)
1519
1520    uas = var0*var3 - var1*var2
1521    vas = var0*var2 + var1*var3
1522
1523    dnamesvar = ['Time','south_north','west_east']
1524    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1525
1526    return uas, vas, dnamesvar, dvnamesvar
1527
1528def compute_WRFta(t, p, dimns, dimvns):
1529    """ Function to compute WRF air temperature
1530      t= orginal WRF temperature
1531      p= original WRF pressure (P + PB)
1532      formula:
1533          temp = theta*(p/p0)**(R/Cp)
1534
1535    """
1536    fname = 'compute_WRFta'
1537
1538    ta = (t+300.)*(p/fdef.module_definitions.p0ref)**(fdef.module_definitions.rcp)
1539
1540    dnamesvar = ['Time','bottom_top','south_north','west_east']
1541    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1542
1543    return ta, dnamesvar, dvnamesvar
1544
1545def compute_WRFtd(t, p, qv, dimns, dimvns):
1546    """ Function to compute WRF dew-point air temperature
1547      t= orginal WRF temperature
1548      p= original WRF pressure (P + PB)
1549      formula:
1550          temp = theta*(p/p0)**(R/Cp)
1551
1552    """
1553    fname = 'compute_WRFtd'
1554
1555    tk = (t+300.)*(p/fdef.module_definitions.p0ref)**(fdef.module_definitions.rcp)
1556    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
1557    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
1558
1559    rh = qv/data2
1560               
1561    pa = rh * data1
1562    td = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
1563
1564    dnamesvar = ['Time','bottom_top','south_north','west_east']
1565    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1566
1567    return td, dnamesvar, dvnamesvar
1568
1569def compute_WRFwd(u, v, sina, cosa, dimns, dimvns):
1570    """ Function to compute the wind direction
1571      u= W-E wind direction [ms-1]
1572      v= N-S wind direction [ms-1]
1573      sina= original WRF local sinus of map rotation
1574      cosa= original WRF local cosinus of map rotation
1575    """
1576    fname = 'compute_WRFwd'
1577    var0 = u
1578    var1 = v
1579    var2 = sina
1580    var3 = cosa
1581
1582    # un-staggering variables
1583    unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
1584    ua = np.zeros(tuple(unstgdims), dtype=np.float)
1585    va = np.zeros(tuple(unstgdims), dtype=np.float)
1586    unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
1587    unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
1588    unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
1589    unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
1590
1591    for iz in range(var0.shape[1]):
1592        ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
1593        va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
1594
1595    theta = np.arctan2(va,ua)
1596    theta = np.where(theta < 0., theta + 2.*np.pi, theta)
1597
1598    wd = 360.*theta/(2.*np.pi)
1599
1600    dnamesvar = ['Time','bottom_top','south_north','west_east']
1601    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1602
1603    return wd
1604
1605####### Variables (as they arrive without dimensions staff)
1606
1607def var_hur(p, t, q):
1608    """ Function to compute relative humidity following 'August - Roche - Magnus' formula
1609      [t]= temperature (assuming [[t],z,y,x] in [K])
1610      [p] = pressure field (assuming in [Pa])
1611      [q] = mixing ratio in [kgkg-1]
1612      ref.: M. G. Lawrence, BAMS, 2005, 225
1613    >>> print var_rh(101300., 290., 3.)
1614    0.250250256174
1615    """
1616    fname = 'var_hur'
1617
1618    # Enthalpy of vaporization [Jkg-1]
1619    L = 2.453*10.**6
1620    # Gas constant for water vapor [JK-1Kg-1]
1621    Rv = 461.5 
1622
1623    # Computing saturated pressure
1624    #es = 10.*0.6112*np.exp(17.67*(t-273.16)/(t-29.65))
1625    #es = 6.11*np.exp(L/Rv*(1.-273./t)/273.)
1626
1627    # August - Roche - Magnus formula
1628    es = 6.1094*np.exp(17.625*(t-273.15)/((t-273.15)+243.04))
1629
1630    # Saturated mixing ratio [g/kg]
1631    ws = 0.622*es/(0.01*p-es)
1632
1633    hur = q/(ws*1000.)
1634
1635    #print 'q:', q[5,5,5,5], 't:', t[5,5,5,5], 'p:', p[5,5,5,5]
1636    #print 'es:', es[5,5,5,5], 'ws:', ws[5,5,5,5], 'hur:', hur[5,5,5,5]
1637
1638    return hur
1639
1640def var_hur_Uhus(p, t, hus):
1641    """ Function to compute relative humidity following 'August - Roche - Magnus' formula using hus
1642      [t]= temperature (assuming [[t],z,y,x] in [K])
1643      [p] = pressure field (assuming in [Pa])
1644      [hus] = specific humidty [1]
1645      ref.: M. G. Lawrence, BAMS, 2005, 225
1646    >>> print var_rh(101300., 290., 3.)
1647    0.250250256174
1648    """
1649    fname = 'var_hur_Uhus'
1650
1651    # Enthalpy of vaporization [Jkg-1]
1652    L = 2.453*10.**6
1653    # Gas constant for water vapor [JK-1Kg-1]
1654    Rv = 461.5 
1655
1656    # Computing saturated pressure
1657    #es = 10.*0.6112*np.exp(17.67*(t-273.16)/(t-29.65))
1658    #es = 6.11*np.exp(L/Rv*(1.-273./t)/273.)
1659
1660    # August - Roche - Magnus formula
1661    es = 6.1094*np.exp(17.625*(t-273.15)/((t-273.15)+243.04))
1662
1663    # Saturated mixing ratio [g/kg]
1664    ws = 0.622*es/(0.01*p-es)
1665
1666    # Mixing ratio
1667    q = hus/(1.+hus)
1668
1669    hur = q/ws
1670
1671    #print 'q:', q[5,5,5,5], 't:', t[5,5,5,5], 'p:', p[5,5,5,5]
1672    #print 'es:', es[5,5,5,5], 'ws:', ws[5,5,5,5], 'hur:', hur[5,5,5,5]
1673
1674    return hur
1675
1676def var_td(t, p, qv):
1677    """ Function to compute dew-point air temperature from temperature and pressure values
1678      t= temperature [K]
1679      p= pressure (Pa)
1680      formula:
1681          temp = theta*(p/p0)**(R/Cp)
1682
1683    """
1684    fname = 'compute_td'
1685
1686    tk = (t)*(p/fdef.module_definitions.p0ref)**(fdef.module_definitions.rcp)
1687    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
1688    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
1689
1690    rh = qv/data2
1691               
1692    pa = rh * data1
1693    td = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
1694
1695    return td
1696
1697def var_td_Uhus(t, p, hus):
1698    """ Function to compute dew-point air temperature from temperature and pressure values using hus
1699      t= temperature [K]
1700      hus= specific humidity [1]
1701      formula:
1702          temp = theta*(p/p0)**(R/Cp)
1703
1704    """
1705    fname = 'compute_td'
1706
1707    tk = (t)*(p/fdef.module_definitions.p0ref)**(fdef.module_definitions.rcp)
1708    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
1709    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
1710
1711    qv = hus/(1.+hus)
1712
1713    rh = qv/data2
1714               
1715    pa = rh * data1
1716    td = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
1717
1718    return td
1719
1720def var_wd(u, v):
1721    """ Function to compute the wind direction
1722      [u]= W-E wind direction [ms-1, knot, ...]
1723      [v]= N-S wind direction [ms-1, knot, ...]
1724    """
1725    fname = 'var_wd'
1726
1727    theta = np.arctan2(v,u)
1728    theta = np.where(theta < 0., theta + 2.*np.pi, theta)
1729
1730    wd = 360.*theta/(2.*np.pi)
1731
1732    return wd
1733
1734def var_ws(u, v):
1735    """ Function to compute the wind speed
1736      [u]= W-E wind direction [ms-1, knot, ...]
1737      [v]= N-S wind direction [ms-1, knot, ...]
1738    """
1739    fname = 'var_ws'
1740
1741    ws = np.sqrt(u*u + v*v)
1742
1743    return ws
1744
1745class C_diagnostic(object):
1746    """ Class to compute generic variables
1747      Cdiag: name of the diagnostic to compute
1748      ncobj: netcdf object with data
1749      sfcvars: dictionary with CF equivalencies of surface variables inside file
1750      vars3D: dictionary with CF equivalencies of 3D variables inside file
1751      dictdims: dictionary with CF equivalencies of dimensions inside file
1752        self.values = Values of the diagnostic
1753        self.dims = Dimensions of the diagnostic
1754        self.units = units of the diagnostic
1755        self.incvars = list of variables from the input netCDF object
1756    """ 
1757    def __init__(self, Cdiag, ncobj, sfcvars, vars3D, dictdims):
1758        fname = 'C_diagnostic'
1759        self.values = None
1760        self.dims = None
1761        self.incvars = ncobj.variables
1762        self.units = None
1763
1764        if Cdiag == 'hur':
1765            """ Computing relative humidity
1766            """
1767            vn = 'hur'
1768            CF3Dvars = ['ta', 'plev', 'q']
1769            for v3D in CF3Dvars:
1770                if not vars3D.has_key(v3D):
1771                    print gen.errormsg
1772                    print '  ' + fname + ": missing variable '" + v3D +              \
1773                      "' attribution to compute '" + vn + "' !!"
1774                    print '  Equivalence of 3D variables provided _______'
1775                    gen.printing_dictionary(vars3D)
1776                    quit(-1)
1777                if not self.incvars.has_key(vars3D[v3D]):
1778                    print gen.errormsg
1779                    print '  ' + fname + ": missing variable '" + vars3D[v3D] +      \
1780                      "' in input file to compute '" + vn + "' !!"
1781                    print '  available variables:', self.incvars.keys()
1782                    print '  looking for variables _______' 
1783                    gen.printing_dictionary(vars3D)
1784                    quit(-1)
1785
1786            ta = ncobj.variables[vars3D['ta']][:]
1787            p = ncobj.variables[vars3D['plev']][:]
1788            q = ncobj.variables[vars3D['q']][:]
1789
1790            if len(ta.shape) != len(p.shape):
1791                p = gen.fill_Narray(p, ta*0., filldim=[0,2,3])
1792
1793            self.values = var_hur(p, ta, q)
1794            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1795              dictdims['lon']]
1796            self.units = '1'
1797
1798        if Cdiag == 'hur_Uhus':
1799            """ Computing relative humidity using hus
1800            """
1801            vn = 'hur'
1802            CF3Dvars = ['ta', 'plev', 'hus']
1803            for v3D in CF3Dvars:
1804                if not vars3D.has_key(v3D):
1805                    print gen.errormsg
1806                    print '  ' + fname + ": missing variable '" + v3D +              \
1807                      "' attribution to compute '" + vn + "' !!"
1808                    print '  Equivalence of 3D variables provided _______'
1809                    gen.printing_dictionary(vars3D)
1810                    quit(-1)
1811                if not self.incvars.has_key(vars3D[v3D]):
1812                    print gen.errormsg
1813                    print '  ' + fname + ": missing variable '" + vars3D[v3D] +      \
1814                      "' in input file to compute '" + vn + "' !!"
1815                    print '  available variables:', self.incvars.keys()
1816                    print '  looking for variables _______' 
1817                    gen.printing_dictionary(vars3D)
1818                    quit(-1)
1819
1820            ta = ncobj.variables[vars3D['ta']][:]
1821            p = ncobj.variables[vars3D['plev']][:]
1822            hus = ncobj.variables[vars3D['hus']][:]
1823
1824            if len(ta.shape) != len(p.shape):
1825                p = gen.fill_Narray(p, ta*0., filldim=[0,2,3])
1826
1827            self.values = var_hur_Uhus(p, ta, hus)
1828            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1829              dictdims['lon']]
1830            self.units = '1'
1831
1832        elif Cdiag == 'td':
1833            """ Computing dew-point temperature
1834            """
1835            vn = 'td'
1836            CF3Dvars = ['ta', 'plev', 'q']
1837            for v3D in CF3Dvars:
1838                if not vars3D.has_key(v3D):
1839                    print gen.errormsg
1840                    print '  ' + fname + ": missing variable '" + v3D +              \
1841                      "' attribution to compute '" + vn + "' !!"
1842                    print '  Equivalence of 3D variables provided _______'
1843                    gen.printing_dictionary(vars3D)
1844                    quit(-1)
1845                if not self.incvars.has_key(vars3D[v3D]):
1846                    print gen.errormsg
1847                    print '  ' + fname + ": missing variable '" + vars3D[v3D] +      \
1848                      "' in input file to compute '" + vn + "' !!"
1849                    print '  available variables:', self.incvars.keys()
1850                    print '  looking for variables _______' 
1851                    gen.printing_dictionary(vars3D)
1852                    quit(-1)
1853
1854            ta = ncobj.variables[vars3D['ta']][:]
1855            p = ncobj.variables[vars3D['plev']][:]
1856            q = ncobj.variables[vars3D['q']][:]
1857
1858            if len(ta.shape) != len(p.shape):
1859                p = gen.fill_Narray(p, ta*0., filldim=[0,2,3])
1860
1861            self.values = var_td(ta, p, q)
1862            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1863              dictdims['lon']]
1864            self.units = 'K'
1865
1866        elif Cdiag == 'td_Uhus':
1867            """ Computing dew-point temperature
1868            """
1869            vn = 'td'
1870            CF3Dvars = ['ta', 'plev', 'hus']
1871            for v3D in CF3Dvars:
1872                if not vars3D.has_key(v3D):
1873                    print gen.errormsg
1874                    print '  ' + fname + ": missing variable '" + v3D +              \
1875                      "' attribution to compute '" + vn + "' !!"
1876                    print '  Equivalence of 3D variables provided _______'
1877                    gen.printing_dictionary(vars3D)
1878                    quit(-1)
1879                if not self.incvars.has_key(vars3D[v3D]):
1880                    print gen.errormsg
1881                    print '  ' + fname + ": missing variable '" + vars3D[v3D] +      \
1882                      "' in input file to compute '" + vn + "' !!"
1883                    print '  available variables:', self.incvars.keys()
1884                    print '  looking for variables _______' 
1885                    gen.printing_dictionary(vars3D)
1886                    quit(-1)
1887
1888            ta = ncobj.variables[vars3D['ta']][:]
1889            p = ncobj.variables[vars3D['plev']][:]
1890            hus = ncobj.variables[vars3D['hus']][:]
1891
1892            if len(ta.shape) != len(p.shape):
1893                p = gen.fill_Narray(p, ta*0., filldim=[0,2,3])
1894
1895            self.values = var_td_Uhus(ta, p, hus)
1896            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1897              dictdims['lon']]
1898            self.units = 'K'
1899
1900        elif Cdiag == 'wd':
1901            """ Computing wind direction
1902            """
1903            vn = 'wd'
1904            CF3Dvars = ['ua', 'va']
1905            for v3D in CF3Dvars:
1906                if not vars3D.has_key(v3D):
1907                    print gen.errormsg
1908                    print '  ' + fname + ": missing variable '" + v3D +              \
1909                      "self.' attribution to compute '" + vn + "' !!"
1910                    print '  Equivalence of 3D variables provided _______'
1911                    gen.printing_dictionary(vars3D)
1912                    quit(-1)
1913                if not self.incvars.has_key(vars3D[v3D]):
1914                    print gen.errormsg
1915                    print '  ' + fname + ": missing variable '" + vars3D[v3D] +      \
1916                      "' in input file to compute '" + vn + "' !!"
1917                    print '  available variables:', self.incvars.keys()
1918                    print '  looking for variables _______' 
1919                    gen.printing_dictionary(vars3D)
1920                    quit(-1)
1921   
1922            ua = ncobj.variables[vars3D['ua']][:]
1923            va = ncobj.variables[vars3D['va']][:]
1924   
1925            self.values = var_wd(ua, va)
1926            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1927              dictdims['lon']]
1928            self.units = 'degree'
1929
1930        elif Cdiag == 'ws':
1931            """ Computing wind speed
1932            """
1933            vn = 'ws'
1934            CF3Dvars = ['ua', 'va']
1935            for v3D in CF3Dvars:
1936                if not vars3D.has_key(v3D):
1937                    print gen.errormsg
1938                    print '  ' + fname + ": missing variable '" + v3D +              \
1939                      "' attribution to compute '" + vn + "' !!"
1940                    print '  Equivalence of 3D variables provided _______'
1941                    gen.printing_dictionary(vars3D)
1942                    quit(-1)
1943                if not self.incvars.has_key(vars3D[v3D]):
1944                    print gen.errormsg
1945                    print '  ' + fname + ": missing variable '" + vars3D[v3D] +      \
1946                      "' in input file to compute '" + vn + "' !!"
1947                    print '  available variables:', self.incvars.keys()
1948                    print '  looking for variables _______' 
1949                    gen.printing_dictionary(vars3D)
1950                    quit(-1)
1951
1952            ua = ncobj.variables[vars3D['ua']][:]
1953            va = ncobj.variables[vars3D['va']][:]
1954
1955            self.values = var_ws(ua, va)
1956            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
1957              dictdims['lon']]
1958            self.units = ncobj.variables[vars3D['ua']].units
1959
1960        else:
1961            print gen.errormsg
1962            print '  ' + fname + ": variable '" + Wdiag + "' not ready !!"
1963            print '  available ones:', Cavailablediags
1964            quit(-1)
1965
1966class W_diagnostic(object):
1967    """ Class to compute WRF diagnostics variables
1968      Wdiag: name of the diagnostic to compute
1969      ncobj: netcdf object with data
1970      sfcvars: dictionary with CF equivalencies of surface variables inside file
1971      vars3D: dictionary with CF equivalencies of 3D variables inside file
1972      indims: list of dimensions inside file
1973      invardims: list of dimension-variables inside file
1974      dictdims: dictionary with CF equivalencies of dimensions inside file
1975        self.values = Values of the diagnostic
1976        self.dims = Dimensions of the diagnostic
1977        self.units = units of the diagnostic
1978        self.incvars = list of variables from the input netCDF object
1979    """   
1980    def __init__(self, Wdiag, ncobj, sfcvars, vars3D, indims, invardims, dictdims):
1981        fname = 'W_diagnostic'
1982
1983        self.values = None
1984        self.dims = None
1985        self.incvars = ncobj.variables
1986        self.units = None
1987
1988        if Wdiag == 'hur':
1989            """ Computing relative humidity
1990            """
1991            vn = 'hur'
1992            CF3Dvars = ['ta', 'hus']
1993            for v3D in CF3Dvars:
1994                if not vars3D.has_key(v3D):
1995                    print gen.errormsg
1996                    print '  ' + fname + ": missing variable '" + v3D +              \
1997                      "' attribution to compute '" + vn + "' !!"
1998                    print '  Equivalence of 3D variables provided _______'
1999                    gen.printing_dictionary(vars3D)
2000                    quit(-1)
2001
2002            ta = ncobj.variables['T'][:]
2003            p = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
2004            hur = ncobj.variables['QVAPOR'][:]
2005   
2006            vals, dims, vdims = compute_WRFhur(ta, p, hur, indims, invardims)
2007            self.values = vals
2008            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
2009              dictdims['lon']]
2010            self.units = '1'
2011
2012        elif Wdiag == 'p':
2013            """ Computing air pressure
2014            """
2015            vn = 'p'
2016
2017            self.values = ncobj.variables['PB'][:] + ncobj.variables['P'][:]
2018            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
2019              dictdims['lon']]
2020            self.units = ncobj.variables['PB'].units
2021
2022        elif Wdiag == 'ta':
2023            """ Computing air temperature
2024            """
2025            vn = 'ta'
2026            CF3Dvars = ['ta']
2027            for v3D in CF3Dvars:
2028                if not vars3D.has_key(v3D):
2029                    print gen.errormsg
2030                    print '  ' + fname + ": missing variable '" + v3D +              \
2031                      "' attribution to compute '" + vn + "' !!"
2032                    print '  Equivalence of 3D variables provided _______'
2033                    gen.printing_dictionary(vars3D)
2034                    quit(-1)
2035
2036            ta = ncobj.variables['T'][:]
2037            p = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
2038   
2039            vals, dims, vdims = compute_WRFta(ta, p, indims, invardims)
2040            self.values = vals
2041            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
2042              dictdims['lon']]
2043            self.units = 'K'
2044
2045        elif Wdiag == 'td':
2046            """ Computing dew-point temperature
2047            """
2048            vn = 'td'
2049            CF3Dvars = ['ta', 'hus']
2050            for v3D in CF3Dvars:
2051                if not vars3D.has_key(v3D):
2052                    print gen.errormsg
2053                    print '  ' + fname + ": missing variable '" + v3D +              \
2054                      "' attribution to compute '" + vn + "' !!"
2055                    print '  Equivalence of 3D variables provided _______'
2056                    gen.printing_dictionary(vars3D)
2057                    quit(-1)
2058
2059            ta = ncobj.variables['T'][:]
2060            p = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
2061            hus = ncobj.variables['QVAPOR'][:]
2062   
2063            vals, dims, vdims = compute_WRFtd(ta, p, hus, indims, invardims)
2064            self.values = vals
2065            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
2066              dictdims['lon']]
2067            self.units = 'K'
2068   
2069        elif Wdiag == 'ua':
2070            """ Computing x-wind
2071            """
2072            vn = 'ua'
2073            CF3Dvars = ['ua', 'va']
2074            for v3D in CF3Dvars:
2075                if not vars3D.has_key(v3D):
2076                    print gen.errormsg
2077                    print '  ' + fname + ": missing variable '" + v3D +              \
2078                      "' attribution to compute '" + vn + "' !!"
2079                    print '  Equivalence of 3D variables provided _______'
2080                    gen.printing_dictionary(vars3D)
2081                    quit(-1)
2082
2083            ua = ncobj.variables['U'][:]
2084            va = ncobj.variables['V'][:]
2085            sina = ncobj.variables['SINALPHA'][:]
2086            cosa = ncobj.variables['COSALPHA'][:]
2087   
2088            vals, dims, vdims = compute_WRFua(ua, va, sina, cosa, indims, invardims)
2089            self.values = vals
2090            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
2091              dictdims['lon']]
2092            self.units = ncobj.variables['U'].units
2093   
2094        elif Wdiag == 'uas':
2095            """ Computing 10m x-wind
2096            """
2097            vn = 'uas'
2098            CFsfcvars = ['uas', 'vas']
2099            for vsf in CFsfcvars:
2100                if not sfcvars.has_key(vsf):
2101                    print gen.errormsg
2102                    print '  ' + fname + ": missing variable '" + vsf +              \
2103                      "' attribution to compute '" + vn + "' !!"
2104                    print '  Equivalence of sfc variables provided _______'
2105                    gen.printing_dictionary(sfcvars)
2106                    quit(-1)
2107   
2108            uas = ncobj.variables['U10'][:]
2109            vas = ncobj.variables['V10'][:]
2110            sina = ncobj.variables['SINALPHA'][:]
2111            cosa = ncobj.variables['COSALPHA'][:]
2112   
2113            vals,dims,vdims = compute_WRFuas(uas, vas, sina, cosa, indims, invardims)
2114            self.values = vals
2115            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
2116              dictdims['lon']]
2117            self.units = ncobj.variables['U10'].units
2118   
2119        elif Wdiag == 'va':
2120            """ Computing y-wind
2121            """
2122            vn = 'ua'
2123            CF3Dvars = ['ua', 'va']
2124            for v3D in CF3Dvars:
2125                if not vars3D.has_key(v3D):
2126                    print gen.errormsg
2127                    print '  ' + fname + ": missing variable '" + v3D +              \
2128                      "' attribution to compute '" + vn + "' !!"
2129                    print '  Equivalence of 3D variables provided _______'
2130                    gen.printing_dictionary(vars3D)
2131                    quit(-1)
2132   
2133            ua = ncobj.variables['U'][:]
2134            va = ncobj.variables['V'][:]
2135            sina = ncobj.variables['SINALPHA'][:]
2136            cosa = ncobj.variables['COSALPHA'][:]
2137   
2138            vals, dims, vdims = compute_WRFva(ua, va, sina, cosa, indims, invardims)
2139            self.values = vals
2140            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
2141              dictdims['lon']]
2142            self.units = ncobj.variables['U'].units
2143   
2144        elif Wdiag == 'vas':
2145            """ Computing 10m y-wind
2146            """
2147            vn = 'uas'
2148            CFsfcvars = ['uas', 'vas']
2149            for vsf in CFsfcvars:
2150                if not sfcvars.has_key(vsf):
2151                    print gen.errormsg
2152                    print '  ' + fname + ": missing variable '" + vsf +              \
2153                      "' attribution to compute '" + vn + "' !!"
2154                    print '  Equivalence of sfc variables provided _______'
2155                    gen.printing_dictionary(sfcvars)
2156                    quit(-1)
2157   
2158            uas = ncobj.variables['U10'][:]
2159            vas = ncobj.variables['V10'][:]
2160            sina = ncobj.variables['SINALPHA'][:]
2161            cosa = ncobj.variables['COSALPHA'][:]
2162   
2163            vals,dims,vdims = compute_WRFvas(uas, vas, sina, cosa, indims, invardims)
2164            self.values = vals
2165            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
2166              dictdims['lon']]
2167            self.units = ncobj.variables['U10'].units
2168
2169        elif Wdiag == 'wd':
2170            """ Computing wind direction
2171            """
2172            vn = 'wd'
2173            CF3Dvars = ['ua', 'va']
2174            for v3D in CF3Dvars:
2175                if not vars3D.has_key(v3D):
2176                    print gen.errormsg
2177                    print '  ' + fname + ": missing variable '" + v3D +              \
2178                      "' attribution to compute '" + vn + "' !!"
2179                    print '  Equivalence of 3D variables provided _______'
2180                    gen.printing_dictionary(vars3D)
2181                    quit(-1)
2182
2183            ua = ncobj.variables['U10'][:]
2184            va = ncobj.variables['V10'][:]
2185            sina = ncobj.variables['SINALPHA'][:]
2186            cosa = ncobj.variables['COSALPHA'][:]
2187   
2188            vals, dims, vdims = compute_WRFwd(ua, va, sina, cosa, indims, invardims)
2189            self.values = vals
2190            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
2191              dictdims['lon']]
2192            self.units = 'degree'
2193
2194        elif Wdiag == 'ws':
2195            """ Computing wind speed
2196            """
2197            vn = 'ws'
2198            CF3Dvars = ['ua', 'va']
2199            for v3D in CF3Dvars:
2200                if not vars3D.has_key(v3D):
2201                    print gen.errormsg
2202                    print '  ' + fname + ": missing variable '" + v3D +              \
2203                      "' attribution to compute '" + vn + "' !!"
2204                    print '  Equivalence of 3D variables provided _______'
2205                    gen.printing_dictionary(vars3D)
2206                    quit(-1)
2207   
2208            ua = ncobj.variables['U10'][:]
2209            va = ncobj.variables['V10'][:]
2210   
2211            self.values = var_ws(ua, va)
2212            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
2213              dictdims['lon']]
2214            self.units = ncobj.variables['U10'].units
2215
2216        elif Wdiag == 'zg':
2217            """ Computing geopotential
2218            """
2219            vn = 'zg'
2220
2221            self.values = ncobj.variables['PHB'][:] + ncobj.variables['PH'][:]
2222            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
2223              dictdims['lon']]
2224            self.units = ncobj.variables['PHB'].units
2225
2226        else:
2227            print gen.errormsg
2228            print '  ' + fname + ": variable '" + Wdiag + "' not ready !!"
2229            print '  available ones:', Wavailablediags
2230            quit(-1)
Note: See TracBrowser for help on using the repository browser.