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

Last change on this file since 1682 was 1680, checked in by lfita, 8 years ago

Fix typo
File size: 33.3 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# compute_wds: Function to compute the wind direction
19# compute_wss: Function to compute the wind speed
20# compute_WRFta: Function to compute WRF air temperature
21# compute_WRFtd: Function to compute WRF dew-point air temperature
22# compute_WRFuava: Function to compute geographical rotated WRF 3D winds
23# compute_WRFuasvas: Fucntion to compute geographical rotated WRF 2-meter winds
24# derivate_centered: Function to compute the centered derivate of a given field
25# def Forcompute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ via Fortran subroutine
26# Forcompute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ via a Fortran module
27
28# Others just providing variable values
29# var_cllmh: Fcuntion to compute cllmh on a 1D column
30# var_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ using 1D vertical column values
31# var_mslp: Fcuntion to compute mean sea-level pressure
32# var_virtualTemp: This function returns virtual temperature in K,
33# var_WRFtime: Function to copmute CFtimes from WRFtime variable
34# rotational_z: z-component of the rotatinoal of horizontal vectorial field
35# turbulence_var: Function to compute the Taylor's decomposition turbulence term from a a given variable
36
37import numpy as np
38from netCDF4 import Dataset as NetCDFFile
39import os
40import re
41import nc_var_tools as ncvar
42import generic_tools as gen
43import datetime as dtime
44import module_ForDiag as fdin
45import module_ForDef as fdef
46
47main = 'diag_tools.py'
48errormsg = 'ERROR -- error -- ERROR -- error'
49warnmsg = 'WARNING -- warning -- WARNING -- warning'
50
51# Constants
52grav = fdef.module_definitions.grav
53
54# Gneral information
55##
56def reduce_spaces(string):
57    """ Function to give words of a line of text removing any extra space
58    """
59    values = string.replace('\n','').split(' ')
60    vals = []
61    for val in values:
62         if len(val) > 0:
63             vals.append(val)
64
65    return vals
66
67def variable_combo(varn,combofile):
68    """ Function to provide variables combination from a given variable name
69      varn= name of the variable
70      combofile= ASCII file with the combination of variables
71        [varn] [combo]
72          [combo]: '@' separated list of variables to use to generate [varn]
73            [WRFdt] to get WRF time-step (from general attributes)
74    >>> variable_combo('WRFprls','/home/lluis/PY/diagnostics.inf')
75    deaccum@RAINNC@XTIME@prnc
76    """
77    fname = 'variable_combo'
78
79    if varn == 'h':
80        print fname + '_____________________________________________________________'
81        print variable_combo.__doc__
82        quit()
83
84    if not os.path.isfile(combofile):
85        print errormsg
86        print '  ' + fname + ": file with combinations '" + combofile +              \
87          "' does not exist!!"
88        quit(-1)
89
90    objf = open(combofile, 'r')
91
92    found = False
93    for line in objf:
94        linevals = reduce_spaces(line)
95        varnf = linevals[0]
96        combo = linevals[1].replace('\n','')
97        if varn == varnf: 
98            found = True
99            break
100
101    if not found:
102        print errormsg
103        print '  ' + fname + ": variable '" + varn + "' not found in '" + combofile +\
104          "' !!"
105        combo='ERROR'
106
107    objf.close()
108
109    return combo
110
111# Mathematical operators
112##
113def compute_accum(varv, dimns, dimvns):
114    """ Function to compute the accumulation of a variable
115    compute_accum(varv, dimnames, dimvns)
116      [varv]= values to accum (assuming [t,])
117      [dimns]= list of the name of the dimensions of the [varv]
118      [dimvns]= list of the name of the variables with the values of the
119        dimensions of [varv]
120    """
121    fname = 'compute_accum'
122
123    deacdims = dimns[:]
124    deacvdims = dimvns[:]
125
126    slicei = []
127    slicee = []
128
129    Ndims = len(varv.shape)
130    for iid in range(0,Ndims):
131        slicei.append(slice(0,varv.shape[iid]))
132        slicee.append(slice(0,varv.shape[iid]))
133
134    slicee[0] = np.arange(varv.shape[0])
135    slicei[0] = np.arange(varv.shape[0])
136    slicei[0][1:varv.shape[0]] = np.arange(varv.shape[0]-1)
137
138    vari = varv[tuple(slicei)]
139    vare = varv[tuple(slicee)]
140
141    ac = vari*0.
142    for it in range(1,varv.shape[0]):
143        ac[it,] = ac[it-1,] + vare[it,]
144
145    return ac, deacdims, deacvdims
146
147def compute_deaccum(varv, dimns, dimvns):
148    """ Function to compute the deaccumulation of a variable
149    compute_deaccum(varv, dimnames, dimvns)
150      [varv]= values to deaccum (assuming [t,])
151      [dimns]= list of the name of the dimensions of the [varv]
152      [dimvns]= list of the name of the variables with the values of the
153        dimensions of [varv]
154    """
155    fname = 'compute_deaccum'
156
157    deacdims = dimns[:]
158    deacvdims = dimvns[:]
159
160    slicei = []
161    slicee = []
162
163    Ndims = len(varv.shape)
164    for iid in range(0,Ndims):
165        slicei.append(slice(0,varv.shape[iid]))
166        slicee.append(slice(0,varv.shape[iid]))
167
168    slicee[0] = np.arange(varv.shape[0])
169    slicei[0] = np.arange(varv.shape[0])
170    slicei[0][1:varv.shape[0]] = np.arange(varv.shape[0]-1)
171
172    vari = varv[tuple(slicei)]
173    vare = varv[tuple(slicee)]
174
175    deac = vare - vari
176
177    return deac, deacdims, deacvdims
178
179def derivate_centered(var,dim,dimv):
180    """ Function to compute the centered derivate of a given field
181      centered derivate(n) = (var(n-1) + var(n+1))/(2*dn).
182    [var]= variable
183    [dim]= which dimension to compute the derivate
184    [dimv]= dimension values (can be of different dimension of [var])
185    >>> derivate_centered(np.arange(16).reshape(4,4)*1.,1,1.)
186    [[  0.   1.   2.   0.]
187     [  0.   5.   6.   0.]
188     [  0.   9.  10.   0.]
189     [  0.  13.  14.   0.]]
190    """
191
192    fname = 'derivate_centered'
193   
194    vark = var.dtype
195
196    if hasattr(dimv, "__len__"):
197# Assuming that the last dimensions of var [..., N, M] are the same of dimv [N, M]
198        if len(var.shape) != len(dimv.shape):
199            dimvals = np.zeros((var.shape), dtype=vark)
200            if len(var.shape) - len(dimv.shape) == 1:
201                for iz in range(var.shape[0]):
202                    dimvals[iz,] = dimv
203            elif len(var.shape) - len(dimv.shape) == 2:
204                for it in range(var.shape[0]):
205                    for iz in range(var.shape[1]):
206                        dimvals[it,iz,] = dimv
207            else:
208                print errormsg
209                print '  ' + fname + ': dimension difference between variable',      \
210                  var.shape,'and variable with dimension values',dimv.shape,         \
211                  ' not ready !!!'
212                quit(-1)
213        else:
214            dimvals = dimv
215    else:
216# dimension values are identical everywhere!
217# from: http://stackoverflow.com/questions/16807011/python-how-to-identify-if-a-variable-is-an-array-or-a-scalar   
218        dimvals = np.ones((var.shape), dtype=vark)*dimv
219
220    derivate = np.zeros((var.shape), dtype=vark)
221    if dim > len(var.shape) - 1:
222        print errormsg
223        print '  ' + fname + ': dimension',dim,' too big for given variable of ' +   \
224          'shape:', var.shape,'!!!'
225        quit(-1)
226
227    slicebef = []
228    sliceaft = []
229    sliceder = []
230
231    for id in range(len(var.shape)):
232        if id == dim:
233            slicebef.append(slice(0,var.shape[id]-2))
234            sliceaft.append(slice(2,var.shape[id]))
235            sliceder.append(slice(1,var.shape[id]-1))
236        else:
237            slicebef.append(slice(0,var.shape[id]))
238            sliceaft.append(slice(0,var.shape[id]))
239            sliceder.append(slice(0,var.shape[id]))
240
241    if hasattr(dimv, "__len__"):
242        derivate[tuple(sliceder)] = (var[tuple(slicebef)] + var[tuple(sliceaft)])/   \
243          ((dimvals[tuple(sliceaft)] - dimvals[tuple(slicebef)]))
244        print (dimvals[tuple(sliceaft)] - dimvals[tuple(slicebef)])
245    else:
246        derivate[tuple(sliceder)] = (var[tuple(slicebef)] + var[tuple(sliceaft)])/   \
247          (2.*dimv)
248
249#    print 'before________'
250#    print var[tuple(slicebef)]
251
252#    print 'after________'
253#    print var[tuple(sliceaft)]
254
255    return derivate
256
257def rotational_z(Vx,Vy,pos):
258    """ z-component of the rotatinoal of horizontal vectorial field
259    \/ x (Vx,Vy,Vz) = \/xVy - \/yVx
260    [Vx]= Variable component x
261    [Vy]=  Variable component y
262    [pos]= poisition of the grid points
263    >>> rotational_z(np.arange(16).reshape(4,4)*1., np.arange(16).reshape(4,4)*1., 1.)
264    [[  0.   1.   2.   0.]
265     [ -4.   0.   0.  -7.]
266     [ -8.   0.   0. -11.]
267     [  0.  13.  14.   0.]]
268    """
269
270    fname =  'rotational_z'
271
272    ndims = len(Vx.shape)
273    rot1 = derivate_centered(Vy,ndims-1,pos)
274    rot2 = derivate_centered(Vx,ndims-2,pos)
275
276    rot = rot1 - rot2
277
278    return rot
279
280# Diagnostics
281##
282
283def var_clt(cfra):
284    """ Function to compute the total cloud fraction following 'newmicro.F90' from
285      LMDZ using 1D vertical column values
286      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
287    """
288    ZEPSEC=1.0E-12
289
290    fname = 'var_clt'
291
292    zclear = 1.
293    zcloud = 0.
294
295    dz = cfra.shape[0]
296    for iz in range(dz):
297        zclear =zclear*(1.-np.max([cfra[iz],zcloud]))/(1.-np.min([zcloud,1.-ZEPSEC]))
298        clt = 1. - zclear
299        zcloud = cfra[iz]
300
301    return clt
302
303def compute_clt(cldfra, dimns, dimvns):
304    """ Function to compute the total cloud fraction following 'newmicro.F90' from
305      LMDZ
306    compute_clt(cldfra, dimnames)
307      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
308      [dimns]= list of the name of the dimensions of [cldfra]
309      [dimvns]= list of the name of the variables with the values of the
310        dimensions of [cldfra]
311    """
312    fname = 'compute_clt'
313
314    cltdims = dimns[:]
315    cltvdims = dimvns[:]
316
317    if len(cldfra.shape) == 4:
318        clt = np.zeros((cldfra.shape[0],cldfra.shape[2],cldfra.shape[3]),            \
319          dtype=np.float)
320        dx = cldfra.shape[3]
321        dy = cldfra.shape[2]
322        dz = cldfra.shape[1]
323        dt = cldfra.shape[0]
324        cltdims.pop(1)
325        cltvdims.pop(1)
326
327        for it in range(dt):
328            for ix in range(dx):
329                for iy in range(dy):
330                    zclear = 1.
331                    zcloud = 0.
332                    gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
333                    clt[it,iy,ix] = var_clt(cldfra[it,:,iy,ix])
334
335    else:
336        clt = np.zeros((cldfra.shape[1],cldfra.shape[2]), dtype=np.float)
337        dx = cldfra.shape[2]
338        dy = cldfra.shape[1]
339        dy = cldfra.shape[0]
340        cltdims.pop(0)
341        cltvdims.pop(0)
342        for ix in range(dx):
343            for iy in range(dy):
344                zclear = 1.
345                zcloud = 0.
346                gen.percendone(ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
347                clt[iy,ix] = var_clt(cldfra[:,iy,ix])
348
349    return clt, cltdims, cltvdims
350
351def Forcompute_clt(cldfra, dimns, dimvns):
352    """ Function to compute the total cloud fraction following 'newmicro.F90' from
353      LMDZ via a Fortran module
354    compute_clt(cldfra, dimnames)
355      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
356      [dimns]= list of the name of the dimensions of [cldfra]
357      [dimvns]= list of the name of the variables with the values of the
358        dimensions of [cldfra]
359    """
360    fname = 'Forcompute_clt'
361
362    cltdims = dimns[:]
363    cltvdims = dimvns[:]
364
365
366    if len(cldfra.shape) == 4:
367        clt = np.zeros((cldfra.shape[0],cldfra.shape[2],cldfra.shape[3]),            \
368          dtype=np.float)
369        dx = cldfra.shape[3]
370        dy = cldfra.shape[2]
371        dz = cldfra.shape[1]
372        dt = cldfra.shape[0]
373        cltdims.pop(1)
374        cltvdims.pop(1)
375
376        clt = fdin.module_fordiagnostics.compute_clt4d2(cldfra[:])
377
378    else:
379        clt = np.zeros((cldfra.shape[1],cldfra.shape[2]), dtype=np.float)
380        dx = cldfra.shape[2]
381        dy = cldfra.shape[1]
382        dy = cldfra.shape[0]
383        cltdims.pop(0)
384        cltvdims.pop(0)
385
386        clt = fdin.module_fordiagnostics.compute_clt3d1(cldfra[:])
387
388    return clt, cltdims, cltvdims
389
390def var_cllmh(cfra, p):
391    """ Fcuntion to compute cllmh on a 1D column
392    """
393
394    fname = 'var_cllmh'
395
396    ZEPSEC =1.0E-12
397    prmhc = 440.*100.
398    prmlc = 680.*100.
399
400    zclearl = 1.
401    zcloudl = 0.
402    zclearm = 1.
403    zcloudm = 0.
404    zclearh = 1.
405    zcloudh = 0.
406
407    dvz = cfra.shape[0]
408
409    cllmh = np.ones((3), dtype=np.float)
410
411    for iz in range(dvz):
412        if p[iz] < prmhc:
413            cllmh[2] = cllmh[2]*(1.-np.max([cfra[iz], zcloudh]))/(1.-                \
414              np.min([zcloudh,1.-ZEPSEC]))
415            zcloudh = cfra[iz]
416        elif p[iz] >= prmhc and p[iz] < prmlc:
417            cllmh[1] = cllmh[1]*(1.-np.max([cfra[iz], zcloudm]))/(1.-                \
418              np.min([zcloudm,1.-ZEPSEC]))
419            zcloudm = cfra[iz]
420        elif p[iz] >= prmlc:
421            cllmh[0] = cllmh[0]*(1.-np.max([cfra[iz], zcloudl]))/(1.-                \
422              np.min([zcloudl,1.-ZEPSEC]))
423            zcloudl = cfra[iz]
424
425    cllmh = 1.- cllmh
426
427    return cllmh
428
429def Forcompute_cllmh(cldfra, pres, dimns, dimvns):
430    """ Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ via Fortran subroutine
431    compute_clt(cldfra, pres, dimns, dimvns)
432      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
433      [pres] = pressure field
434      [dimns]= list of the name of the dimensions of [cldfra]
435      [dimvns]= list of the name of the variables with the values of the
436        dimensions of [cldfra]
437    """
438    fname = 'Forcompute_cllmh'
439
440    cllmhdims = dimns[:]
441    cllmhvdims = dimvns[:]
442
443    if len(cldfra.shape) == 4:
444        dx = cldfra.shape[3]
445        dy = cldfra.shape[2]
446        dz = cldfra.shape[1]
447        dt = cldfra.shape[0]
448        cllmhdims.pop(1)
449        cllmhvdims.pop(1)
450
451        cllmh = fdin.module_fordiagnostics.compute_cllmh4d2(cldfra[:], pres[:])
452       
453    else:
454        dx = cldfra.shape[2]
455        dy = cldfra.shape[1]
456        dz = cldfra.shape[0]
457        cllmhdims.pop(0)
458        cllmhvdims.pop(0)
459
460        cllmh = fdin.module_fordiagnostics.compute_cllmh3d1(cldfra[:], pres[:])
461
462    return cllmh, cllmhdims, cllmhvdims
463
464def compute_cllmh(cldfra, pres, dimns, dimvns):
465    """ Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ
466    compute_clt(cldfra, pres, dimns, dimvns)
467      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
468      [pres] = pressure field
469      [dimns]= list of the name of the dimensions of [cldfra]
470      [dimvns]= list of the name of the variables with the values of the
471        dimensions of [cldfra]
472    """
473    fname = 'compute_cllmh'
474
475    cllmhdims = dimns[:]
476    cllmhvdims = dimvns[:]
477
478    if len(cldfra.shape) == 4:
479        dx = cldfra.shape[3]
480        dy = cldfra.shape[2]
481        dz = cldfra.shape[1]
482        dt = cldfra.shape[0]
483        cllmhdims.pop(1)
484        cllmhvdims.pop(1)
485
486        cllmh = np.ones(tuple([3, dt, dy, dx]), dtype=np.float)
487
488        for it in range(dt):
489            for ix in range(dx):
490                for iy in range(dy):
491                    gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
492                    cllmh[:,it,iy,ix] = var_cllmh(cldfra[it,:,iy,ix], pres[it,:,iy,ix])
493       
494    else:
495        dx = cldfra.shape[2]
496        dy = cldfra.shape[1]
497        dz = cldfra.shape[0]
498        cllmhdims.pop(0)
499        cllmhvdims.pop(0)
500
501        cllmh = np.ones(tuple([3, dy, dx]), dtype=np.float)
502
503        for ix in range(dx):
504            for iy in range(dy):
505                gen.percendone(ix*dy + iy,dx*dy, 5, 'diagnosted')
506                cllmh[:,iy,ix] = var_cllmh(cldfra[:,iy,ix], pres[:,iy,ix])
507
508    return cllmh, cllmhdims, cllmhvdims
509
510def compute_clivi(dens, qtot, dimns, dimvns):
511    """ Function to compute cloud-ice water path (clivi)
512      [dens] = density [in kgkg-1] (assuming [t],z,y,x)
513      [qtot] = added mixing ratio of all cloud-ice species in [kgkg-1] (assuming [t],z,y,x)
514      [dimns]= list of the name of the dimensions of [q]
515      [dimvns]= list of the name of the variables with the values of the
516        dimensions of [q]
517    """
518    fname = 'compute_clivi'
519
520    clividims = dimns[:]
521    clivivdims = dimvns[:]
522
523    if len(qtot.shape) == 4:
524        clividims.pop(1)
525        clivivdims.pop(1)
526    else:
527        clividims.pop(0)
528        clivivdims.pop(0)
529
530    data1 = dens*qtot
531    clivi = np.sum(data1, axis=1)
532
533    return clivi, clividims, clivivdims
534
535
536def compute_clwvl(dens, qtot, dimns, dimvns):
537    """ Function to compute condensed water path (clwvl)
538      [dens] = density [in kgkg-1] (assuming [t],z,y,x)
539      [qtot] = added mixing ratio of all cloud-water species in [kgkg-1] (assuming [t],z,y,x)
540      [dimns]= list of the name of the dimensions of [q]
541      [dimvns]= list of the name of the variables with the values of the
542        dimensions of [q]
543    """
544    fname = 'compute_clwvl'
545
546    clwvldims = dimns[:]
547    clwvlvdims = dimvns[:]
548
549    if len(qtot.shape) == 4:
550        clwvldims.pop(1)
551        clwvlvdims.pop(1)
552    else:
553        clwvldims.pop(0)
554        clwvlvdims.pop(0)
555
556    data1 = dens*qtot
557    clwvl = np.sum(data1, axis=1)
558
559    return clwvl, clwvldims, clwvlvdims
560
561def var_virtualTemp (temp,rmix):
562    """ This function returns virtual temperature in K,
563      temp: temperature [K]
564      rmix: mixing ratio in [kgkg-1]
565    """
566
567    fname = 'var_virtualTemp'
568
569    virtual=temp*(0.622+rmix)/(0.622*(1.+rmix))
570
571    return virtual
572
573
574def var_mslp(pres, psfc, ter, tk, qv):
575    """ Function to compute mslp on a 1D column
576    """
577
578    fname = 'var_mslp'
579
580    N = 1.0
581    expon=287.04*.0065/9.81
582    pref = 40000.
583
584# First find where about 400 hPa is located
585    dz=len(pres) 
586
587    kref = -1
588    pinc = pres[0] - pres[dz-1]
589
590    if pinc < 0.:
591        for iz in range(1,dz):
592            if pres[iz-1] >= pref and pres[iz] < pref: 
593                kref = iz
594                break
595    else:
596        for iz in range(dz-1):
597            if pres[iz] >= pref and pres[iz+1] < pref: 
598                kref = iz
599                break
600
601    if kref == -1:
602        print errormsg
603        print '  ' + fname + ': no reference pressure:',pref,'found!!'
604        print '    values:',pres[:]
605        quit(-1)
606
607    mslp = 0.
608
609# We are below both the ground and the lowest data level.
610
611# First, find the model level that is closest to a "target" pressure
612# level, where the "target" pressure is delta-p less that the local
613# value of a horizontally smoothed surface pressure field.  We use
614# delta-p = 150 hPa here. A standard lapse rate temperature profile
615# passing through the temperature at this model level will be used
616# to define the temperature profile below ground.  This is similar
617# to the Benjamin and Miller (1990) method, using 
618# 700 hPa everywhere for the "target" pressure.
619
620# ptarget = psfc - 15000.
621    ptarget = 70000.
622    dpmin=1.e4
623    kupper = 0
624    if pinc > 0.:
625        for iz in range(dz-1,0,-1):
626            kupper = iz
627            dp=np.abs( pres[iz] - ptarget )
628            if dp < dpmin: exit
629            dpmin = np.min([dpmin, dp])
630    else:
631        for iz in range(dz):
632            kupper = iz
633            dp=np.abs( pres[iz] - ptarget )
634            if dp < dpmin: exit
635            dpmin = np.min([dpmin, dp])
636
637    pbot=np.max([pres[0], psfc])
638#    zbot=0.
639
640#    tbotextrap=tk(i,j,kupper,itt)*(pbot/pres_field(i,j,kupper,itt))**expon
641#    tvbotextrap=virtual(tbotextrap,qv(i,j,1,itt))
642
643#    data_out(i,j,itt,1) = (zbot+tvbotextrap/.0065*(1.-(interp_levels(1)/pbot)**expon))
644    tbotextrap = tk[kupper]*(psfc/ptarget)**expon
645    tvbotextrap = var_virtualTemp(tbotextrap, qv[kupper])
646    mslp = psfc*( (tvbotextrap+0.0065*ter)/tvbotextrap)**(1./expon)
647
648    return mslp
649
650def compute_mslp(pressure, psurface, terrain, temperature, qvapor, dimns, dimvns):
651    """ Function to compute mslp: mean sea level pressure following p_interp.F90 from WRF
652    var_mslp(pres, ter, tk, qv, dimns, dimvns)
653      [pressure]= pressure field [Pa] (assuming [[t],z,y,x])
654      [psurface]= surface pressure field [Pa]
655      [terrain]= topography [m]
656      [temperature]= temperature [K]
657      [qvapor]= water vapour mixing ratio [kgkg-1]
658      [dimns]= list of the name of the dimensions of [cldfra]
659      [dimvns]= list of the name of the variables with the values of the
660        dimensions of [pres]
661    """
662
663    fname = 'compute_mslp'
664
665    mslpdims = list(dimns[:])
666    mslpvdims = list(dimvns[:])
667
668    if len(pressure.shape) == 4:
669        mslpdims.pop(1)
670        mslpvdims.pop(1)
671    else:
672        mslpdims.pop(0)
673        mslpvdims.pop(0)
674
675    if len(pressure.shape) == 4:
676        dx = pressure.shape[3]
677        dy = pressure.shape[2]
678        dz = pressure.shape[1]
679        dt = pressure.shape[0]
680
681        mslpv = np.zeros(tuple([dt, dy, dx]), dtype=np.float)
682
683# Terrain... to 2D !
684        terval = np.zeros(tuple([dy, dx]), dtype=np.float)
685        if len(terrain.shape) == 3:
686            terval = terrain[0,:,:]
687        else:
688            terval = terrain
689
690        for ix in range(dx):
691            for iy in range(dy):
692                if terval[iy,ix] > 0.:
693                    for it in range(dt):
694                        mslpv[it,iy,ix] = var_mslp(pressure[it,:,iy,ix],             \
695                          psurface[it,iy,ix], terval[iy,ix], temperature[it,:,iy,ix],\
696                          qvapor[it,:,iy,ix])
697
698                        gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
699                else:
700                    mslpv[:,iy,ix] = psurface[:,iy,ix]
701
702    else:
703        dx = pressure.shape[2]
704        dy = pressure.shape[1]
705        dz = pressure.shape[0]
706
707        mslpv = np.zeros(tuple([dy, dx]), dtype=np.float)
708
709# Terrain... to 2D !
710        terval = np.zeros(tuple([dy, dx]), dtype=np.float)
711        if len(terrain.shape) == 3:
712            terval = terrain[0,:,:]
713        else:
714            terval = terrain
715
716        for ix in range(dx):
717            for iy in range(dy):
718                gen.percendone(ix*dy + iy,dx*dy, 5, 'diagnosted')
719                if terval[iy,ix] > 0.:
720                    mslpv[iy,ix] = var_mslp(pressure[:,iy,ix], psurface[iy,ix],          \
721                      terval[iy,ix], temperature[:,iy,ix], qvapor[:,iy,ix])
722                else:
723                    mslpv[iy,ix] = psfc[iy,ix]
724
725    return mslpv, mslpdims, mslpvdims
726
727def compute_OMEGAw(omega, p, t, dimns, dimvns):
728    """ Function to transform OMEGA [Pas-1] to velocities [ms-1]
729      tacking: https://www.ncl.ucar.edu/Document/Functions/Contributed/omega_to_w.shtml
730      [omega] = vertical velocity [in ms-1] (assuming [t],z,y,x)
731      [p] = pressure in [Pa] (assuming [t],z,y,x)
732      [t] = temperature in [K] (assuming [t],z,y,x)
733      [dimns]= list of the name of the dimensions of [q]
734      [dimvns]= list of the name of the variables with the values of the
735        dimensions of [q]
736    """
737    fname = 'compute_OMEGAw'
738
739    rgas = 287.058            # J/(kg-K) => m2/(s2 K)
740    g    = 9.80665            # m/s2
741
742    wdims = dimns[:]
743    wvdims = dimvns[:]
744
745    rho  = p/(rgas*t)         # density => kg/m3
746    w    = -omega/(rho*g)     
747
748    return w, wdims, wvdims
749
750def compute_prw(dens, q, dimns, dimvns):
751    """ Function to compute water vapour path (prw)
752      [dens] = density [in kgkg-1] (assuming [t],z,y,x)
753      [q] = mixing ratio in [kgkg-1] (assuming [t],z,y,x)
754      [dimns]= list of the name of the dimensions of [q]
755      [dimvns]= list of the name of the variables with the values of the
756        dimensions of [q]
757    """
758    fname = 'compute_prw'
759
760    prwdims = dimns[:]
761    prwvdims = dimvns[:]
762
763    if len(q.shape) == 4:
764        prwdims.pop(1)
765        prwvdims.pop(1)
766    else:
767        prwdims.pop(0)
768        prwvdims.pop(0)
769
770    data1 = dens*q
771    prw = np.sum(data1, axis=1)
772
773    return prw, prwdims, prwvdims
774
775def compute_rh(p, t, q, dimns, dimvns):
776    """ Function to compute relative humidity following 'Tetens' equation (T,P) ...'
777      [t]= temperature (assuming [[t],z,y,x] in [K])
778      [p] = pressure field (assuming in [hPa])
779      [q] = mixing ratio in [kgkg-1]
780      [dimns]= list of the name of the dimensions of [t]
781      [dimvns]= list of the name of the variables with the values of the
782        dimensions of [t]
783    """
784    fname = 'compute_rh'
785
786    rhdims = dimns[:]
787    rhvdims = dimvns[:]
788
789    data1 = 10.*0.6112*np.exp(17.67*(t-273.16)/(t-29.65))
790    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
791
792    rh = q/data2
793
794    return rh, rhdims, rhvdims
795
796def compute_td(p, temp, qv, dimns, dimvns):
797    """ Function to compute the dew point temperature
798      [p]= pressure [Pa]
799      [temp]= temperature [C]
800      [qv]= mixing ratio [kgkg-1]
801      [dimns]= list of the name of the dimensions of [p]
802      [dimvns]= list of the name of the variables with the values of the
803        dimensions of [p]
804    """
805    fname = 'compute_td'
806
807#    print '    ' + fname + ': computing dew-point temperature from TS as t and Tetens...'
808# tacking from: http://en.wikipedia.org/wiki/Dew_point
809    tk = temp
810    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
811    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
812
813    rh = qv/data2
814               
815    pa = rh * data1
816    td = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
817
818    tddims = dimns[:]
819    tdvdims = dimvns[:]
820
821    return td, tddims, tdvdims
822
823def var_WRFtime(timewrfv, refdate='19491201000000', tunitsval='minutes'):
824    """ Function to copmute CFtimes from WRFtime variable
825      refdate= [YYYYMMDDMIHHSS] format of reference date
826      tunitsval= CF time units
827      timewrfv= matrix string values of WRF 'Times' variable
828    """
829    fname = 'var_WRFtime'
830
831    yrref=refdate[0:4]
832    monref=refdate[4:6]
833    dayref=refdate[6:8]
834    horref=refdate[8:10]
835    minref=refdate[10:12]
836    secref=refdate[12:14]
837
838    refdateS = yrref + '-' + monref + '-' + dayref + ' ' + horref + ':' + minref +   \
839      ':' + secref
840
841    dt = timewrfv.shape[0]
842    WRFtime = np.zeros((dt), dtype=np.float)
843
844    for it in range(dt):
845        wrfdates = gen.datetimeStr_conversion(timewrfv[it,:],'WRFdatetime', 'matYmdHMS')
846        WRFtime[it] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
847
848    tunits = tunitsval + ' since ' + refdateS
849
850    return WRFtime, tunits
851
852def turbulence_var(varv, dimvn, dimn):
853    """ Function to compute the Taylor's decomposition turbulence term from a a given variable
854      x*=<x^2>_t-(<X>_t)^2
855    turbulence_var(varv,dimn)
856      varv= values of the variable
857      dimvn= names of the dimension of the variable
858      dimn= names of the dimensions (as a dictionary with 'X', 'Y', 'Z', 'T')
859    >>> turbulence_var(np.arange((27)).reshape(3,3,3),['time','y','x'],{'T':'time', 'Y':'y', 'X':'x'})
860    [[ 54.  54.  54.]
861     [ 54.  54.  54.]
862     [ 54.  54.  54.]]
863    """
864    fname = 'turbulence_varv'
865
866    timedimid = dimvn.index(dimn['T'])
867
868    varv2 = varv*varv
869
870    vartmean = np.mean(varv, axis=timedimid)
871    var2tmean = np.mean(varv2, axis=timedimid)
872
873    varvturb = var2tmean - (vartmean*vartmean)
874
875    return varvturb
876
877def compute_turbulence(v, dimns, dimvns):
878    """ Function to compute the rubulence term of the Taylor's decomposition ...'
879      x*=<x^2>_t-(<X>_t)^2
880      [v]= variable (assuming [[t],z,y,x])
881      [dimns]= list of the name of the dimensions of [v]
882      [dimvns]= list of the name of the variables with the values of the
883        dimensions of [v]
884    """
885    fname = 'compute_turbulence'
886
887    turbdims = dimns[:]
888    turbvdims = dimvns[:]
889
890    turbdims.pop(0)
891    turbvdims.pop(0)
892
893    v2 = v*v
894
895    vartmean = np.mean(v, axis=0)
896    var2tmean = np.mean(v2, axis=0)
897
898    turb = var2tmean - (vartmean*vartmean)
899
900    return turb, turbdims, turbvdims
901
902def compute_wds(u, v, dimns, dimvns):
903    """ Function to compute the wind direction
904      [u]= W-E wind direction [ms-1, knot, ...]
905      [v]= N-S wind direction [ms-1, knot, ...]
906      [dimns]= list of the name of the dimensions of [u]
907      [dimvns]= list of the name of the variables with the values of the
908        dimensions of [u]
909    """
910    fname = 'compute_wds'
911
912#    print '    ' + fname + ': computing wind direction as ATAN2(v,u) ...'
913    theta = np.arctan2(v,u)
914    theta = np.where(theta < 0., theta + 2.*np.pi, theta)
915
916    wds = 360.*theta/(2.*np.pi)
917
918    wdsdims = dimns[:]
919    wdsvdims = dimvns[:]
920
921    return wds, wdsdims, wdsvdims
922
923def compute_wss(u, v, dimns, dimvns):
924    """ Function to compute the wind speed
925      [u]= W-E wind direction [ms-1, knot, ...]
926      [v]= N-S wind direction [ms-1, knot, ...]
927      [dimns]= list of the name of the dimensions of [u]
928      [dimvns]= list of the name of the variables with the values of the
929        dimensions of [u]
930    """
931    fname = 'compute_wss'
932
933#    print '    ' + fname + ': computing wind speed as SQRT(v**2 + u**2) ...'
934    wss = np.sqrt(u*u + v*v)
935
936    wssdims = dimns[:]
937    wssvdims = dimvns[:]
938
939    return wss, wssdims, wssvdims
940
941def timeunits_seconds(dtu):
942    """ Function to transform a time units to seconds
943    timeunits_seconds(timeuv)
944      [dtu]= time units value to transform in seconds
945    """
946    fname='timunits_seconds'
947
948    if dtu == 'years':
949        times = 365.*24.*3600.
950    elif dtu == 'weeks':
951        times = 7.*24.*3600.
952    elif dtu == 'days':
953        times = 24.*3600.
954    elif dtu == 'hours':
955        times = 3600.
956    elif dtu == 'minutes':
957        times = 60.
958    elif dtu == 'seconds':
959        times = 1.
960    elif dtu == 'miliseconds':
961        times = 1./1000.
962    else:
963        print errormsg
964        print '  ' + fname  + ": time units '" + dtu + "' not ready !!"
965        quit(-1)
966
967    return times
968
969def compute_WRFuava(u, v, sina, cosa, dimns, dimvns):
970    """ Function to compute geographical rotated WRF 3D winds
971      u= orginal WRF x-wind
972      v= orginal WRF y-wind
973      sina= original WRF local sinus of map rotation
974      cosa= original WRF local cosinus of map rotation
975      formula:
976        ua = u*cosa-va*sina
977        va = u*sina+va*cosa
978    """
979    fname = 'compute_WRFuava'
980
981    var0 = u
982    var1 = v
983    var2 = sina
984    var3 = cosa
985
986    # un-staggering variables
987    unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
988    ua = np.zeros(tuple(unstgdims), dtype=np.float)
989    va = np.zeros(tuple(unstgdims), dtype=np.float)
990    unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
991    unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
992    unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
993    unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
994
995    for iz in range(var0.shape[1]):
996        ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
997        va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3
998
999    dnamesvar = ['Time','bottom_top','south_north','west_east']
1000    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1001
1002    return ua, va, dnamesvar, dvnamesvar
1003
1004def compute_WRFuasvas(u10, v10, sina, cosa, dimns, dimvns):
1005    """ Function to compute geographical rotated WRF 2-meter winds
1006      u10= orginal WRF 10m x-wind
1007      v10= orginal WRF 10m y-wind
1008      sina= original WRF local sinus of map rotation
1009      cosa= original WRF local cosinus of map rotation
1010      formula:
1011        uas = u10*cosa-va10*sina
1012        vas = u10*sina+va10*cosa
1013    """
1014    fname = 'compute_WRFuasvas'
1015
1016    var0 = u10
1017    var1 = v10
1018    var2 = sina
1019    var3 = cosa
1020
1021    uas = np.zeros(var0.shape, dtype=np.float)
1022    vas = np.zeros(var0.shape, dtype=np.float)
1023
1024    uas = var0*var3 - var1*var2
1025    vas = var0*var2 + var1*var3
1026
1027    dnamesvar = ['Time','south_north','west_east']
1028    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1029
1030    return uas, vas, dnamesvar, dvnamesvar
1031
1032def compute_WRFta(t, p, dimns, dimvns):
1033    """ Function to compute WRF air temperature
1034      t= orginal WRF temperature
1035      p= original WRF pressure (P + PB)
1036      formula:
1037          temp = theta*(p/p0)**(R/Cp)
1038
1039    """
1040    fname = 'compute_WRFta'
1041
1042    ta = (t+300.)*(p/fdef.module_definitions.p0ref)**(fdef.module_definitions.rcp)
1043
1044    dnamesvar = ['Time','bottom_top','south_north','west_east']
1045    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1046
1047    return ta, dnamesvar, dvnamesvar
1048
1049def compute_WRFtd(t, p, qv, dimns, dimvns):
1050    """ Function to compute WRF dew-point air temperature
1051      t= orginal WRF temperature
1052      p= original WRF pressure (P + PB)
1053      formula:
1054          temp = theta*(p/p0)**(R/Cp)
1055
1056    """
1057    fname = 'compute_WRFtd'
1058
1059    tk = (t+300.)*(p/fdef.module_definitions.p0ref)**(fdef.module_definitions.rcp)
1060    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
1061    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
1062
1063    rh = qv/data2
1064               
1065    pa = rh * data1
1066    td = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
1067
1068    dnamesvar = ['Time','bottom_top','south_north','west_east']
1069    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)
1070
1071    return td, dnamesvar, dvnamesvar
Note: See TracBrowser for help on using the repository browser.