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

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

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