Changeset 1675 in lmdz_wrf


Ignore:
Timestamp:
Dec 4, 2017, 5:41:09 PM (7 years ago)
Author:
lfita
Message:

Advancing diagnostics by creation of an independent module with all the variables called 'diag_tools.py'

Location:
trunk/tools
Files:
1 added
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/diagnostics.py

    r1647 r1675  
    1 # Pthon script to comput diagnostics
     1# Python script to comput diagnostics
    22# L. Fita, LMD. CNR, UPMC-Jussieu, Paris, France
    33# File diagnostics.inf provides the combination of variables to get the desired diagnostic
    … …  
    88## e.g. # diagnostics.py -d 'Time@time,bottom_top@ZNU,south_north@XLAT,west_east@XLONG' -v 'clt|CLDFRA,cllmh|CLDFRA@WRFp,RAINTOT|RAINC@RAINNC@XTIME' -f WRF_LMDZ/NPv31/wrfout_d01_1980-03-01_00:00:00
    99## e.g. # diagnostics.py -f /home/lluis/PY/diagnostics.inf -d variable_combo -v WRFprc
     10
     11# Available general pupose diagnostics (model independent) providing (varv1, varv2, ..., dimns, dimvns)
     12# compute_accum: Function to compute the accumulation of a variable
     13# compute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following
     14#   newmicro.F90 from LMDZ compute_clt(cldfra, pres, dimns, dimvns)
     15# compute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ
     16# compute_clivi: Function to compute cloud-ice water path (clivi)
     17# compute_clwvl: Function to compute condensed water path (clwvl)
     18# compute_deaccum: Function to compute the deaccumulation of a variable
     19# compute_mslp: Function to compute mslp: mean sea level pressure following p_interp.F90 from WRF
     20# compute_OMEGAw: Function to transform OMEGA [Pas-1] to velocities [ms-1]
     21# compute_prw: Function to compute water vapour path (prw)
     22# compute_rh: Function to compute relative humidity following 'Tetens' equation (T,P) ...'
     23# compute_td: Function to compute the dew point temperature
     24# compute_turbulence: Function to compute the rubulence term of the Taylor's decomposition ...'
     25# compute_wds: Function to compute the wind direction
     26# compute_wss: Function to compute the wind speed
     27# compute_WRFuava: Function to compute geographical rotated WRF 3D winds
     28# compute_WRFuasvas: Fucntion to compute geographical rotated WRF 2-meter winds
     29# derivate_centered: Function to compute the centered derivate of a given field
     30# def Forcompute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ via Fortran subroutine
     31# Forcompute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ via a Fortran module
     32
     33# Others just providing variable values
     34# var_cllmh: Fcuntion to compute cllmh on a 1D column
     35# var_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ using 1D vertical column values
     36# var_mslp: Fcuntion to compute mean sea-level pressure
     37# var_virtualTemp: This function returns virtual temperature in K,
     38# var_WRFtime: Function to copmute CFtimes from WRFtime variable
     39# rotational_z: z-component of the rotatinoal of horizontal vectorial field
     40# turbulence_var: Function to compute the Taylor's decomposition turbulence term from a a given variable
    1041
    1142from optparse import OptionParser
    … …  
    1849import datetime as dtime
    1950import module_ForDiag as fdin
     51import diag_tools as diag
    2052
    2153main = 'diagnostics.py'
    … …  
    2658grav = 9.81
    2759
    28 # Gneral information
    29 ##
    30 def reduce_spaces(string):
    31     """ Function to give words of a line of text removing any extra space
    32     """
    33     values = string.replace('\n','').split(' ')
    34     vals = []
    35     for val in values:
    36          if len(val) > 0:
    37              vals.append(val)
    38 
    39     return vals
    40 
    41 def variable_combo(varn,combofile):
    42     """ Function to provide variables combination from a given variable name
    43       varn= name of the variable
    44       combofile= ASCII file with the combination of variables
    45         [varn] [combo]
    46           [combo]: '@' separated list of variables to use to generate [varn]
    47             [WRFdt] to get WRF time-step (from general attributes)
    48     >>> variable_combo('WRFprls','/home/lluis/PY/diagnostics.inf')
    49     deaccum@RAINNC@XTIME@prnc
    50     """
    51     fname = 'variable_combo'
    52 
    53     if varn == 'h':
    54         print fname + '_____________________________________________________________'
    55         print variable_combo.__doc__
    56         quit()
    57 
    58     if not os.path.isfile(combofile):
    59         print errormsg
    60         print '  ' + fname + ": file with combinations '" + combofile +              \
    61           "' does not exist!!"
    62         quit(-1)
    63 
    64     objf = open(combofile, 'r')
    65 
    66     found = False
    67     for line in objf:
    68         linevals = reduce_spaces(line)
    69         varnf = linevals[0]
    70         combo = linevals[1].replace('\n','')
    71         if varn == varnf:
    72             found = True
    73             break
    74 
    75     if not found:
    76         print errormsg
    77         print '  ' + fname + ": variable '" + varn + "' not found in '" + combofile +\
    78           "' !!"
    79         combo='ERROR'
    80 
    81     objf.close()
    82 
    83     return combo
    84 
    85 # Mathematical operators
    86 ##
    87 def compute_accum(varv, dimns, dimvns):
    88     """ Function to compute the accumulation of a variable
    89     compute_accum(varv, dimnames, dimvns)
    90       [varv]= values to accum (assuming [t,])
    91       [dimns]= list of the name of the dimensions of the [varv]
    92       [dimvns]= list of the name of the variables with the values of the
    93         dimensions of [varv]
    94     """
    95     fname = 'compute_accum'
    96 
    97     deacdims = dimns[:]
    98     deacvdims = dimvns[:]
    99 
    100     slicei = []
    101     slicee = []
    102 
    103     Ndims = len(varv.shape)
    104     for iid in range(0,Ndims):
    105         slicei.append(slice(0,varv.shape[iid]))
    106         slicee.append(slice(0,varv.shape[iid]))
    107 
    108     slicee[0] = np.arange(varv.shape[0])
    109     slicei[0] = np.arange(varv.shape[0])
    110     slicei[0][1:varv.shape[0]] = np.arange(varv.shape[0]-1)
    111 
    112     vari = varv[tuple(slicei)]
    113     vare = varv[tuple(slicee)]
    114 
    115     ac = vari*0.
    116     for it in range(1,varv.shape[0]):
    117         ac[it,] = ac[it-1,] + vare[it,]
    118 
    119     return ac, deacdims, deacvdims
    120 
    121 def compute_deaccum(varv, dimns, dimvns):
    122     """ Function to compute the deaccumulation of a variable
    123     compute_deaccum(varv, dimnames, dimvns)
    124       [varv]= values to deaccum (assuming [t,])
    125       [dimns]= list of the name of the dimensions of the [varv]
    126       [dimvns]= list of the name of the variables with the values of the
    127         dimensions of [varv]
    128     """
    129     fname = 'compute_deaccum'
    130 
    131     deacdims = dimns[:]
    132     deacvdims = dimvns[:]
    133 
    134     slicei = []
    135     slicee = []
    136 
    137     Ndims = len(varv.shape)
    138     for iid in range(0,Ndims):
    139         slicei.append(slice(0,varv.shape[iid]))
    140         slicee.append(slice(0,varv.shape[iid]))
    141 
    142     slicee[0] = np.arange(varv.shape[0])
    143     slicei[0] = np.arange(varv.shape[0])
    144     slicei[0][1:varv.shape[0]] = np.arange(varv.shape[0]-1)
    145 
    146     vari = varv[tuple(slicei)]
    147     vare = varv[tuple(slicee)]
    148 
    149     deac = vare - vari
    150 
    151     return deac, deacdims, deacvdims
    152 
    153 def derivate_centered(var,dim,dimv):
    154     """ Function to compute the centered derivate of a given field
    155       centered derivate(n) = (var(n-1) + var(n+1))/(2*dn).
    156     [var]= variable
    157     [dim]= which dimension to compute the derivate
    158     [dimv]= dimension values (can be of different dimension of [var])
    159     >>> derivate_centered(np.arange(16).reshape(4,4)*1.,1,1.)
    160     [[  0.   1.   2.   0.]
    161      [  0.   5.   6.   0.]
    162      [  0.   9.  10.   0.]
    163      [  0.  13.  14.   0.]]
    164     """
    165 
    166     fname = 'derivate_centered'
    167    
    168     vark = var.dtype
    169 
    170     if hasattr(dimv, "__len__"):
    171 # Assuming that the last dimensions of var [..., N, M] are the same of dimv [N, M]
    172         if len(var.shape) != len(dimv.shape):
    173             dimvals = np.zeros((var.shape), dtype=vark)
    174             if len(var.shape) - len(dimv.shape) == 1:
    175                 for iz in range(var.shape[0]):
    176                     dimvals[iz,] = dimv
    177             elif len(var.shape) - len(dimv.shape) == 2:
    178                 for it in range(var.shape[0]):
    179                     for iz in range(var.shape[1]):
    180                         dimvals[it,iz,] = dimv
    181             else:
    182                 print errormsg
    183                 print '  ' + fname + ': dimension difference between variable',      \
    184                   var.shape,'and variable with dimension values',dimv.shape,         \
    185                   ' not ready !!!'
    186                 quit(-1)
    187         else:
    188             dimvals = dimv
    189     else:
    190 # dimension values are identical everywhere!
    191 # from: http://stackoverflow.com/questions/16807011/python-how-to-identify-if-a-variable-is-an-array-or-a-scalar   
    192         dimvals = np.ones((var.shape), dtype=vark)*dimv
    193 
    194     derivate = np.zeros((var.shape), dtype=vark)
    195     if dim > len(var.shape) - 1:
    196         print errormsg
    197         print '  ' + fname + ': dimension',dim,' too big for given variable of ' +   \
    198           'shape:', var.shape,'!!!'
    199         quit(-1)
    200 
    201     slicebef = []
    202     sliceaft = []
    203     sliceder = []
    204 
    205     for id in range(len(var.shape)):
    206         if id == dim:
    207             slicebef.append(slice(0,var.shape[id]-2))
    208             sliceaft.append(slice(2,var.shape[id]))
    209             sliceder.append(slice(1,var.shape[id]-1))
    210         else:
    211             slicebef.append(slice(0,var.shape[id]))
    212             sliceaft.append(slice(0,var.shape[id]))
    213             sliceder.append(slice(0,var.shape[id]))
    214 
    215     if hasattr(dimv, "__len__"):
    216         derivate[tuple(sliceder)] = (var[tuple(slicebef)] + var[tuple(sliceaft)])/   \
    217           ((dimvals[tuple(sliceaft)] - dimvals[tuple(slicebef)]))
    218         print (dimvals[tuple(sliceaft)] - dimvals[tuple(slicebef)])
    219     else:
    220         derivate[tuple(sliceder)] = (var[tuple(slicebef)] + var[tuple(sliceaft)])/   \
    221           (2.*dimv)
    222 
    223 #    print 'before________'
    224 #    print var[tuple(slicebef)]
    225 
    226 #    print 'after________'
    227 #    print var[tuple(sliceaft)]
    228 
    229     return derivate
    230 
    231 def rotational_z(Vx,Vy,pos):
    232     """ z-component of the rotatinoal of horizontal vectorial field
    233     \/ x (Vx,Vy,Vz) = \/xVy - \/yVx
    234     [Vx]= Variable component x
    235     [Vy]=  Variable component y
    236     [pos]= poisition of the grid points
    237     >>> rotational_z(np.arange(16).reshape(4,4)*1., np.arange(16).reshape(4,4)*1., 1.)
    238     [[  0.   1.   2.   0.]
    239      [ -4.   0.   0.  -7.]
    240      [ -8.   0.   0. -11.]
    241      [  0.  13.  14.   0.]]
    242     """
    243 
    244     fname =  'rotational_z'
    245 
    246     ndims = len(Vx.shape)
    247     rot1 = derivate_centered(Vy,ndims-1,pos)
    248     rot2 = derivate_centered(Vx,ndims-2,pos)
    249 
    250     rot = rot1 - rot2
    251 
    252     return rot
    253 
    254 # Diagnostics
    255 ##
    256 
    257 def var_clt(cfra):
    258     """ Function to compute the total cloud fraction following 'newmicro.F90' from
    259       LMDZ using 1D vertical column values
    260       [cldfra]= cloud fraction values (assuming [[t],z,y,x])
    261     """
    262     ZEPSEC=1.0E-12
    263 
    264     fname = 'var_clt'
    265 
    266     zclear = 1.
    267     zcloud = 0.
    268 
    269     dz = cfra.shape[0]
    270     for iz in range(dz):
    271         zclear =zclear*(1.-np.max([cfra[iz],zcloud]))/(1.-np.min([zcloud,1.-ZEPSEC]))
    272         clt = 1. - zclear
    273         zcloud = cfra[iz]
    274 
    275     return clt
    276 
    277 def compute_clt(cldfra, dimns, dimvns):
    278     """ Function to compute the total cloud fraction following 'newmicro.F90' from
    279       LMDZ
    280     compute_clt(cldfra, dimnames)
    281       [cldfra]= cloud fraction values (assuming [[t],z,y,x])
    282       [dimns]= list of the name of the dimensions of [cldfra]
    283       [dimvns]= list of the name of the variables with the values of the
    284         dimensions of [cldfra]
    285     """
    286     fname = 'compute_clt'
    287 
    288     cltdims = dimns[:]
    289     cltvdims = dimvns[:]
    290 
    291     if len(cldfra.shape) == 4:
    292         clt = np.zeros((cldfra.shape[0],cldfra.shape[2],cldfra.shape[3]),            \
    293           dtype=np.float)
    294         dx = cldfra.shape[3]
    295         dy = cldfra.shape[2]
    296         dz = cldfra.shape[1]
    297         dt = cldfra.shape[0]
    298         cltdims.pop(1)
    299         cltvdims.pop(1)
    300 
    301         for it in range(dt):
    302             for ix in range(dx):
    303                 for iy in range(dy):
    304                     zclear = 1.
    305                     zcloud = 0.
    306                     gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
    307                     clt[it,iy,ix] = var_clt(cldfra[it,:,iy,ix])
    308 
    309     else:
    310         clt = np.zeros((cldfra.shape[1],cldfra.shape[2]), dtype=np.float)
    311         dx = cldfra.shape[2]
    312         dy = cldfra.shape[1]
    313         dy = cldfra.shape[0]
    314         cltdims.pop(0)
    315         cltvdims.pop(0)
    316         for ix in range(dx):
    317             for iy in range(dy):
    318                 zclear = 1.
    319                 zcloud = 0.
    320                 gen.percendone(ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
    321                 clt[iy,ix] = var_clt(cldfra[:,iy,ix])
    322 
    323     return clt, cltdims, cltvdims
    324 
    325 def Forcompute_clt(cldfra, dimns, dimvns):
    326     """ Function to compute the total cloud fraction following 'newmicro.F90' from
    327       LMDZ via a Fortran module
    328     compute_clt(cldfra, dimnames)
    329       [cldfra]= cloud fraction values (assuming [[t],z,y,x])
    330       [dimns]= list of the name of the dimensions of [cldfra]
    331       [dimvns]= list of the name of the variables with the values of the
    332         dimensions of [cldfra]
    333     """
    334     fname = 'Forcompute_clt'
    335 
    336     cltdims = dimns[:]
    337     cltvdims = dimvns[:]
    338 
    339 
    340     if len(cldfra.shape) == 4:
    341         clt = np.zeros((cldfra.shape[0],cldfra.shape[2],cldfra.shape[3]),            \
    342           dtype=np.float)
    343         dx = cldfra.shape[3]
    344         dy = cldfra.shape[2]
    345         dz = cldfra.shape[1]
    346         dt = cldfra.shape[0]
    347         cltdims.pop(1)
    348         cltvdims.pop(1)
    349 
    350         clt = fdin.module_fordiagnostics.compute_clt4d2(cldfra[:])
    351 
    352     else:
    353         clt = np.zeros((cldfra.shape[1],cldfra.shape[2]), dtype=np.float)
    354         dx = cldfra.shape[2]
    355         dy = cldfra.shape[1]
    356         dy = cldfra.shape[0]
    357         cltdims.pop(0)
    358         cltvdims.pop(0)
    359 
    360         clt = fdin.module_fordiagnostics.compute_clt3d1(cldfra[:])
    361 
    362     return clt, cltdims, cltvdims
    363 
    364 def var_cllmh(cfra, p):
    365     """ Fcuntion to compute cllmh on a 1D column
    366     """
    367 
    368     fname = 'var_cllmh'
    369 
    370     ZEPSEC =1.0E-12
    371     prmhc = 440.*100.
    372     prmlc = 680.*100.
    373 
    374     zclearl = 1.
    375     zcloudl = 0.
    376     zclearm = 1.
    377     zcloudm = 0.
    378     zclearh = 1.
    379     zcloudh = 0.
    380 
    381     dvz = cfra.shape[0]
    382 
    383     cllmh = np.ones((3), dtype=np.float)
    384 
    385     for iz in range(dvz):
    386         if p[iz] < prmhc:
    387             cllmh[2] = cllmh[2]*(1.-np.max([cfra[iz], zcloudh]))/(1.-                \
    388               np.min([zcloudh,1.-ZEPSEC]))
    389             zcloudh = cfra[iz]
    390         elif p[iz] >= prmhc and p[iz] < prmlc:
    391             cllmh[1] = cllmh[1]*(1.-np.max([cfra[iz], zcloudm]))/(1.-                \
    392               np.min([zcloudm,1.-ZEPSEC]))
    393             zcloudm = cfra[iz]
    394         elif p[iz] >= prmlc:
    395             cllmh[0] = cllmh[0]*(1.-np.max([cfra[iz], zcloudl]))/(1.-                \
    396               np.min([zcloudl,1.-ZEPSEC]))
    397             zcloudl = cfra[iz]
    398 
    399     cllmh = 1.- cllmh
    400 
    401     return cllmh
    402 
    403 def Forcompute_cllmh(cldfra, pres, dimns, dimvns):
    404     """ Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ via Fortran subroutine
    405     compute_clt(cldfra, pres, dimns, dimvns)
    406       [cldfra]= cloud fraction values (assuming [[t],z,y,x])
    407       [pres] = pressure field
    408       [dimns]= list of the name of the dimensions of [cldfra]
    409       [dimvns]= list of the name of the variables with the values of the
    410         dimensions of [cldfra]
    411     """
    412     fname = 'Forcompute_cllmh'
    413 
    414     cllmhdims = dimns[:]
    415     cllmhvdims = dimvns[:]
    416 
    417     if len(cldfra.shape) == 4:
    418         dx = cldfra.shape[3]
    419         dy = cldfra.shape[2]
    420         dz = cldfra.shape[1]
    421         dt = cldfra.shape[0]
    422         cllmhdims.pop(1)
    423         cllmhvdims.pop(1)
    424 
    425         cllmh = fdin.module_fordiagnostics.compute_cllmh4d2(cldfra[:], pres[:])
    426        
    427     else:
    428         dx = cldfra.shape[2]
    429         dy = cldfra.shape[1]
    430         dz = cldfra.shape[0]
    431         cllmhdims.pop(0)
    432         cllmhvdims.pop(0)
    433 
    434         cllmh = fdin.module_fordiagnostics.compute_cllmh3d1(cldfra[:], pres[:])
    435 
    436     return cllmh, cllmhdims, cllmhvdims
    437 
    438 def compute_cllmh(cldfra, pres, dimns, dimvns):
    439     """ Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ
    440     compute_clt(cldfra, pres, dimns, dimvns)
    441       [cldfra]= cloud fraction values (assuming [[t],z,y,x])
    442       [pres] = pressure field
    443       [dimns]= list of the name of the dimensions of [cldfra]
    444       [dimvns]= list of the name of the variables with the values of the
    445         dimensions of [cldfra]
    446     """
    447     fname = 'compute_cllmh'
    448 
    449     cllmhdims = dimns[:]
    450     cllmhvdims = dimvns[:]
    451 
    452     if len(cldfra.shape) == 4:
    453         dx = cldfra.shape[3]
    454         dy = cldfra.shape[2]
    455         dz = cldfra.shape[1]
    456         dt = cldfra.shape[0]
    457         cllmhdims.pop(1)
    458         cllmhvdims.pop(1)
    459 
    460         cllmh = np.ones(tuple([3, dt, dy, dx]), dtype=np.float)
    461 
    462         for it in range(dt):
    463             for ix in range(dx):
    464                 for iy in range(dy):
    465                     gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
    466                     cllmh[:,it,iy,ix] = var_cllmh(cldfra[it,:,iy,ix], pres[it,:,iy,ix])
    467        
    468     else:
    469         dx = cldfra.shape[2]
    470         dy = cldfra.shape[1]
    471         dz = cldfra.shape[0]
    472         cllmhdims.pop(0)
    473         cllmhvdims.pop(0)
    474 
    475         cllmh = np.ones(tuple([3, dy, dx]), dtype=np.float)
    476 
    477         for ix in range(dx):
    478             for iy in range(dy):
    479                 gen.percendone(ix*dy + iy,dx*dy, 5, 'diagnosted')
    480                 cllmh[:,iy,ix] = var_cllmh(cldfra[:,iy,ix], pres[:,iy,ix])
    481 
    482     return cllmh, cllmhdims, cllmhvdims
    483 
    484 def compute_clivi(dens, qtot, dimns, dimvns):
    485     """ Function to compute cloud-ice water path (clivi)
    486       [dens] = density [in kgkg-1] (assuming [t],z,y,x)
    487       [qtot] = added mixing ratio of all cloud-ice species in [kgkg-1] (assuming [t],z,y,x)
    488       [dimns]= list of the name of the dimensions of [q]
    489       [dimvns]= list of the name of the variables with the values of the
    490         dimensions of [q]
    491     """
    492     fname = 'compute_clivi'
    493 
    494     clividims = dimns[:]
    495     clivivdims = dimvns[:]
    496 
    497     if len(qtot.shape) == 4:
    498         clividims.pop(1)
    499         clivivdims.pop(1)
    500     else:
    501         clividims.pop(0)
    502         clivivdims.pop(0)
    503 
    504     data1 = dens*qtot
    505     clivi = np.sum(data1, axis=1)
    506 
    507     return clivi, clividims, clivivdims
    508 
    509 
    510 def compute_clwvl(dens, qtot, dimns, dimvns):
    511     """ Function to compute condensed water path (clwvl)
    512       [dens] = density [in kgkg-1] (assuming [t],z,y,x)
    513       [qtot] = added mixing ratio of all cloud-water 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_clwvl'
    519 
    520     clwvldims = dimns[:]
    521     clwvlvdims = dimvns[:]
    522 
    523     if len(qtot.shape) == 4:
    524         clwvldims.pop(1)
    525         clwvlvdims.pop(1)
    526     else:
    527         clwvldims.pop(0)
    528         clwvlvdims.pop(0)
    529 
    530     data1 = dens*qtot
    531     clwvl = np.sum(data1, axis=1)
    532 
    533     return clwvl, clwvldims, clwvlvdims
    534 
    535 def var_virtualTemp (temp,rmix):
    536     """ This function returns virtual temperature in K,
    537       temp: temperature [K]
    538       rmix: mixing ratio in [kgkg-1]
    539     """
    540 
    541     fname = 'var_virtualTemp'
    542 
    543     virtual=temp*(0.622+rmix)/(0.622*(1.+rmix))
    544 
    545     return virtual
    546 
    547 
    548 def var_mslp(pres, psfc, ter, tk, qv):
    549     """ Function to compute mslp on a 1D column
    550     """
    551 
    552     fname = 'var_mslp'
    553 
    554     N = 1.0
    555     expon=287.04*.0065/9.81
    556     pref = 40000.
    557 
    558 # First find where about 400 hPa is located
    559     dz=len(pres)
    560 
    561     kref = -1
    562     pinc = pres[0] - pres[dz-1]
    563 
    564     if pinc < 0.:
    565         for iz in range(1,dz):
    566             if pres[iz-1] >= pref and pres[iz] < pref:
    567                 kref = iz
    568                 break
    569     else:
    570         for iz in range(dz-1):
    571             if pres[iz] >= pref and pres[iz+1] < pref:
    572                 kref = iz
    573                 break
    574 
    575     if kref == -1:
    576         print errormsg
    577         print '  ' + fname + ': no reference pressure:',pref,'found!!'
    578         print '    values:',pres[:]
    579         quit(-1)
    580 
    581     mslp = 0.
    582 
    583 # We are below both the ground and the lowest data level.
    584 
    585 # First, find the model level that is closest to a "target" pressure
    586 # level, where the "target" pressure is delta-p less that the local
    587 # value of a horizontally smoothed surface pressure field.  We use
    588 # delta-p = 150 hPa here. A standard lapse rate temperature profile
    589 # passing through the temperature at this model level will be used
    590 # to define the temperature profile below ground.  This is similar
    591 # to the Benjamin and Miller (1990) method, using 
    592 # 700 hPa everywhere for the "target" pressure.
    593 
    594 # ptarget = psfc - 15000.
    595     ptarget = 70000.
    596     dpmin=1.e4
    597     kupper = 0
    598     if pinc > 0.:
    599         for iz in range(dz-1,0,-1):
    600             kupper = iz
    601             dp=np.abs( pres[iz] - ptarget )
    602             if dp < dpmin: exit
    603             dpmin = np.min([dpmin, dp])
    604     else:
    605         for iz in range(dz):
    606             kupper = iz
    607             dp=np.abs( pres[iz] - ptarget )
    608             if dp < dpmin: exit
    609             dpmin = np.min([dpmin, dp])
    610 
    611     pbot=np.max([pres[0], psfc])
    612 #    zbot=0.
    613 
    614 #    tbotextrap=tk(i,j,kupper,itt)*(pbot/pres_field(i,j,kupper,itt))**expon
    615 #    tvbotextrap=virtual(tbotextrap,qv(i,j,1,itt))
    616 
    617 #    data_out(i,j,itt,1) = (zbot+tvbotextrap/.0065*(1.-(interp_levels(1)/pbot)**expon))
    618     tbotextrap = tk[kupper]*(psfc/ptarget)**expon
    619     tvbotextrap = var_virtualTemp(tbotextrap, qv[kupper])
    620     mslp = psfc*( (tvbotextrap+0.0065*ter)/tvbotextrap)**(1./expon)
    621 
    622     return mslp
    623 
    624 def compute_mslp(pressure, psurface, terrain, temperature, qvapor, dimns, dimvns):
    625     """ Function to compute mslp: mean sea level pressure following p_interp.F90 from WRF
    626     var_mslp(pres, ter, tk, qv, dimns, dimvns)
    627       [pressure]= pressure field [Pa] (assuming [[t],z,y,x])
    628       [psurface]= surface pressure field [Pa]
    629       [terrain]= topography [m]
    630       [temperature]= temperature [K]
    631       [qvapor]= water vapour mixing ratio [kgkg-1]
    632       [dimns]= list of the name of the dimensions of [cldfra]
    633       [dimvns]= list of the name of the variables with the values of the
    634         dimensions of [pres]
    635     """
    636 
    637     fname = 'compute_mslp'
    638 
    639     mslpdims = list(dimns[:])
    640     mslpvdims = list(dimvns[:])
    641 
    642     if len(pressure.shape) == 4:
    643         mslpdims.pop(1)
    644         mslpvdims.pop(1)
    645     else:
    646         mslpdims.pop(0)
    647         mslpvdims.pop(0)
    648 
    649     if len(pressure.shape) == 4:
    650         dx = pressure.shape[3]
    651         dy = pressure.shape[2]
    652         dz = pressure.shape[1]
    653         dt = pressure.shape[0]
    654 
    655         mslpv = np.zeros(tuple([dt, dy, dx]), dtype=np.float)
    656 
    657 # Terrain... to 2D !
    658         terval = np.zeros(tuple([dy, dx]), dtype=np.float)
    659         if len(terrain.shape) == 3:
    660             terval = terrain[0,:,:]
    661         else:
    662             terval = terrain
    663 
    664         for ix in range(dx):
    665             for iy in range(dy):
    666                 if terval[iy,ix] > 0.:
    667                     for it in range(dt):
    668                         mslpv[it,iy,ix] = var_mslp(pressure[it,:,iy,ix],             \
    669                           psurface[it,iy,ix], terval[iy,ix], temperature[it,:,iy,ix],\
    670                           qvapor[it,:,iy,ix])
    671 
    672                         gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
    673                 else:
    674                     mslpv[:,iy,ix] = psurface[:,iy,ix]
    675 
    676     else:
    677         dx = pressure.shape[2]
    678         dy = pressure.shape[1]
    679         dz = pressure.shape[0]
    680 
    681         mslpv = np.zeros(tuple([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                 gen.percendone(ix*dy + iy,dx*dy, 5, 'diagnosted')
    693                 if terval[iy,ix] > 0.:
    694                     mslpv[iy,ix] = var_mslp(pressure[:,iy,ix], psurface[iy,ix],          \
    695                       terval[iy,ix], temperature[:,iy,ix], qvapor[:,iy,ix])
    696                 else:
    697                     mslpv[iy,ix] = psfc[iy,ix]
    698 
    699     return mslpv, mslpdims, mslpvdims
    700 
    701 def compute_OMEGAw(omega, p, t, dimns, dimvns):
    702     """ Function to transform OMEGA [Pas-1] to velocities [ms-1]
    703       tacking: https://www.ncl.ucar.edu/Document/Functions/Contributed/omega_to_w.shtml
    704       [omega] = vertical velocity [in ms-1] (assuming [t],z,y,x)
    705       [p] = pressure in [Pa] (assuming [t],z,y,x)
    706       [t] = temperature in [K] (assuming [t],z,y,x)
    707       [dimns]= list of the name of the dimensions of [q]
    708       [dimvns]= list of the name of the variables with the values of the
    709         dimensions of [q]
    710     """
    711     fname = 'compute_OMEGAw'
    712 
    713     rgas = 287.058            # J/(kg-K) => m2/(s2 K)
    714     g    = 9.80665            # m/s2
    715 
    716     wdims = dimns[:]
    717     wvdims = dimvns[:]
    718 
    719     rho  = p/(rgas*t)         # density => kg/m3
    720     w    = -omega/(rho*g)     
    721 
    722     return w, wdims, wvdims
    723 
    724 def compute_prw(dens, q, dimns, dimvns):
    725     """ Function to compute water vapour path (prw)
    726       [dens] = density [in kgkg-1] (assuming [t],z,y,x)
    727       [q] = mixing ratio in [kgkg-1] (assuming [t],z,y,x)
    728       [dimns]= list of the name of the dimensions of [q]
    729       [dimvns]= list of the name of the variables with the values of the
    730         dimensions of [q]
    731     """
    732     fname = 'compute_prw'
    733 
    734     prwdims = dimns[:]
    735     prwvdims = dimvns[:]
    736 
    737     if len(q.shape) == 4:
    738         prwdims.pop(1)
    739         prwvdims.pop(1)
    740     else:
    741         prwdims.pop(0)
    742         prwvdims.pop(0)
    743 
    744     data1 = dens*q
    745     prw = np.sum(data1, axis=1)
    746 
    747     return prw, prwdims, prwvdims
    748 
    749 def compute_rh(p, t, q, dimns, dimvns):
    750     """ Function to compute relative humidity following 'Tetens' equation (T,P) ...'
    751       [t]= temperature (assuming [[t],z,y,x] in [K])
    752       [p] = pressure field (assuming in [hPa])
    753       [q] = mixing ratio in [kgkg-1]
    754       [dimns]= list of the name of the dimensions of [t]
    755       [dimvns]= list of the name of the variables with the values of the
    756         dimensions of [t]
    757     """
    758     fname = 'compute_rh'
    759 
    760     rhdims = dimns[:]
    761     rhvdims = dimvns[:]
    762 
    763     data1 = 10.*0.6112*np.exp(17.67*(t-273.16)/(t-29.65))
    764     data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
    765 
    766     rh = q/data2
    767 
    768     return rh, rhdims, rhvdims
    769 
    770 def compute_td(p, temp, qv, dimns, dimvns):
    771     """ Function to compute the dew point temperature
    772       [p]= pressure [Pa]
    773       [temp]= temperature [C]
    774       [qv]= mixing ratio [kgkg-1]
    775       [dimns]= list of the name of the dimensions of [p]
    776       [dimvns]= list of the name of the variables with the values of the
    777         dimensions of [p]
    778     """
    779     fname = 'compute_td'
    780 
    781 #    print '    ' + fname + ': computing dew-point temperature from TS as t and Tetens...'
    782 # tacking from: http://en.wikipedia.org/wiki/Dew_point
    783     tk = temp
    784     data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
    785     data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
    786 
    787     rh = qv/data2
    788                
    789     pa = rh * data1
    790     td = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
    791 
    792     tddims = dimns[:]
    793     tdvdims = dimvns[:]
    794 
    795     return td, tddims, tdvdims
    796 
    797 def turbulence_var(varv, dimvn, dimn):
    798     """ Function to compute the Taylor's decomposition turbulence term from a a given variable
    799       x*=<x^2>_t-(<X>_t)^2
    800     turbulence_var(varv,dimn)
    801       varv= values of the variable
    802       dimvn= names of the dimension of the variable
    803       dimn= names of the dimensions (as a dictionary with 'X', 'Y', 'Z', 'T')
    804     >>> turbulence_var(np.arange((27)).reshape(3,3,3),['time','y','x'],{'T':'time', 'Y':'y', 'X':'x'})
    805     [[ 54.  54.  54.]
    806      [ 54.  54.  54.]
    807      [ 54.  54.  54.]]
    808     """
    809     fname = 'turbulence_varv'
    810 
    811     timedimid = dimvn.index(dimn['T'])
    812 
    813     varv2 = varv*varv
    814 
    815     vartmean = np.mean(varv, axis=timedimid)
    816     var2tmean = np.mean(varv2, axis=timedimid)
    817 
    818     varvturb = var2tmean - (vartmean*vartmean)
    819 
    820     return varvturb
    821 
    822 def compute_turbulence(v, dimns, dimvns):
    823     """ Function to compute the rubulence term of the Taylor's decomposition ...'
    824       x*=<x^2>_t-(<X>_t)^2
    825       [v]= variable (assuming [[t],z,y,x])
    826       [dimns]= list of the name of the dimensions of [v]
    827       [dimvns]= list of the name of the variables with the values of the
    828         dimensions of [v]
    829     """
    830     fname = 'compute_turbulence'
    831 
    832     turbdims = dimns[:]
    833     turbvdims = dimvns[:]
    834 
    835     turbdims.pop(0)
    836     turbvdims.pop(0)
    837 
    838     v2 = v*v
    839 
    840     vartmean = np.mean(v, axis=0)
    841     var2tmean = np.mean(v2, axis=0)
    842 
    843     turb = var2tmean - (vartmean*vartmean)
    844 
    845     return turb, turbdims, turbvdims
    846 
    847 def compute_wds(u, v, dimns, dimvns):
    848     """ Function to compute the wind direction
    849       [u]= W-E wind direction [ms-1, knot, ...]
    850       [v]= N-S wind direction [ms-1, knot, ...]
    851       [dimns]= list of the name of the dimensions of [u]
    852       [dimvns]= list of the name of the variables with the values of the
    853         dimensions of [u]
    854     """
    855     fname = 'compute_wds'
    856 
    857 #    print '    ' + fname + ': computing wind direction as ATAN2(v,u) ...'
    858     theta = np.arctan2(v,u)
    859     theta = np.where(theta < 0., theta + 2.*np.pi, theta)
    860 
    861     wds = 360.*theta/(2.*np.pi)
    862 
    863     wdsdims = dimns[:]
    864     wdsvdims = dimvns[:]
    865 
    866     return wds, wdsdims, wdsvdims
    867 
    868 def compute_wss(u, v, dimns, dimvns):
    869     """ Function to compute the wind speed
    870       [u]= W-E wind direction [ms-1, knot, ...]
    871       [v]= N-S wind direction [ms-1, knot, ...]
    872       [dimns]= list of the name of the dimensions of [u]
    873       [dimvns]= list of the name of the variables with the values of the
    874         dimensions of [u]
    875     """
    876     fname = 'compute_wss'
    877 
    878 #    print '    ' + fname + ': computing wind speed as SQRT(v**2 + u**2) ...'
    879     wss = np.sqrt(u*u + v*v)
    880 
    881     wssdims = dimns[:]
    882     wssvdims = dimvns[:]
    883 
    884     return wss, wssdims, wssvdims
    885 
    886 def timeunits_seconds(dtu):
    887     """ Function to transform a time units to seconds
    888     timeunits_seconds(timeuv)
    889       [dtu]= time units value to transform in seconds
    890     """
    891     fname='timunits_seconds'
    892 
    893     if dtu == 'years':
    894         times = 365.*24.*3600.
    895     elif dtu == 'weeks':
    896         times = 7.*24.*3600.
    897     elif dtu == 'days':
    898         times = 24.*3600.
    899     elif dtu == 'hours':
    900         times = 3600.
    901     elif dtu == 'minutes':
    902         times = 60.
    903     elif dtu == 'seconds':
    904         times = 1.
    905     elif dtu == 'miliseconds':
    906         times = 1./1000.
    907     else:
    908         print errormsg
    909         print '  ' + fname  + ": time units '" + dtu + "' not ready !!"
    910         quit(-1)
    911 
    912     return times
    91360
    91461####### ###### ##### #### ### ## #
    … …  
    1239386        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
    1240387
    1241         diagout, diagoutd, diagoutvd = compute_accum(var0,dnamesvar,dvnamesvar)
     388        diagout, diagoutd, diagoutvd = diag.compute_accum(var0,dnamesvar,dvnamesvar)
    1242389
    1243390        CFvarn = ncvar.variables_values(depvars[0])[0]
    … …  
    1268415                var1 = var01
    1269416
    1270         diagout, diagoutd, diagoutvd = Forcompute_cllmh(var0,var1,dnames,dvnames)
     417        diagout, diagoutd, diagoutvd = diag.Forcompute_cllmh(var0,var1,dnames,dvnames)
    1271418
    1272419        # Removing the nonChecking variable-dimensions from the initial list
    … …  
    1284431           
    1285432        var0 = ncobj.variables[depvars]
    1286         diagout, diagoutd, diagoutvd = Forcompute_clt(var0,dnames,dvnames)
     433        diagout, diagoutd, diagoutvd = diag.Forcompute_clt(var0,dnames,dvnames)
    1287434
    1288435        # Removing the nonChecking variable-dimensions from the initial list
    … …  
    1304451        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
    1305452
    1306         diagout, diagoutd, diagoutvd = compute_deaccum(var0,dnamesvar,dvnamesvar)
     453        diagout, diagoutd, diagoutvd = diag.compute_deaccum(var0,dnamesvar,dvnamesvar)
    1307454
    1308455# Transforming to a flux
    … …  
    1324471        var2 = ncobj.variables[depvars[2]][:]
    1325472
    1326         diagout, diagoutd, diagoutvd = compute_rh(var0,var1,var2,dnames,dvnames)
     473        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnames,dvnames)
    1327474        ncvar.insert_variable(ncobj, 'hur', diagout, diagoutd, diagoutvd, newnc)
    1328475
    … …  
    1337484        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
    1338485
    1339         diagout, diagoutd, diagoutvd = compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
     486        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
    1340487
    1341488        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
    … …  
    1359506        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
    1360507
    1361         diagout, diagoutd, diagoutvd = compute_mslp(var0, var1, var2, var3, var4,    \
     508        diagout, diagoutd, diagoutvd = diag.compute_mslp(var0, var1, var2, var3, var4,    \
    1362509          dnamesvar, dvnamesvar)
    1363510
    … …  
    1380527        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
    1381528
    1382         diagout, diagoutd, diagoutvd = compute_OMEGAw(var0,var1,var2,dnamesvar,dvnamesvar)
     529        diagout, diagoutd, diagoutvd = diag.compute_OMEGAw(var0,var1,var2,dnamesvar,dvnamesvar)
    1383530
    1384531        ncvar.insert_variable(ncobj, 'wa', diagout, diagoutd, diagoutvd, newnc)
    … …  
    1399546        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
    1400547
    1401         diagout, diagoutd, diagoutvd = compute_deaccum(var,dnamesvar,dvnamesvar)
     548        diagout, diagoutd, diagoutvd = diag.compute_deaccum(var,dnamesvar,dvnamesvar)
    1402549
    1403550# Transforming to a flux
    … …  
    1448595        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
    1449596
    1450         diagout, diagoutd, diagoutvd = compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
     597        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
    1451598
    1452599        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
    … …  
    1462609        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
    1463610
    1464         diagout, diagoutd, diagoutvd = compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
     611        diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
    1465612
    1466613        ncvar.insert_variable(ncobj, 'tds', diagout, diagoutd, diagoutvd, newnc)
    … …  
    1477624        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
    1478625
    1479         diagout, diagoutd, diagoutvd = compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
     626        diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)
    1480627
    1481628        ncvar.insert_variable(ncobj, 'tds', diagout, diagoutd, diagoutvd, newnc)
    … …  
    1490637        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
    1491638
    1492         diagout, diagoutd, diagoutvd = compute_wds(var0,var1,dnamesvar,dvnamesvar)
     639        diagout, diagoutd, diagoutvd = diag.compute_wds(var0,var1,dnamesvar,dvnamesvar)
    1493640
    1494641        ncvar.insert_variable(ncobj, 'wds', diagout, diagoutd, diagoutvd, newnc)
    … …  
    1503650        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
    1504651
    1505         diagout, diagoutd, diagoutvd = compute_wss(var0,var1,dnamesvar,dvnamesvar)
     652        diagout, diagoutd, diagoutvd = diag.compute_wss(var0,var1,dnamesvar,dvnamesvar)
    1506653
    1507654        ncvar.insert_variable(ncobj, 'wss', diagout, diagoutd, diagoutvd, newnc)
    … …  
    1515662        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
    1516663
    1517         diagout, diagoutd, diagoutvd = compute_turbulence(var0,dnamesvar,dvnamesvar)
     664        diagout, diagoutd, diagoutvd = diag.compute_turbulence(var0,dnamesvar,dvnamesvar)
    1518665        valsvar = gen.variables_values(depvars)
    1519666
    … …  
    1556703        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
    1557704
    1558         diagout, diagoutd, diagoutvd = compute_clivi(var0, qtotv, dnamesvar,dvnamesvar)
     705        diagout, diagoutd, diagoutvd = diag.compute_clivi(var0, qtotv, dnamesvar,dvnamesvar)
    1559706
    1560707        # Removing the nonChecking variable-dimensions from the initial list
    … …  
    1581728        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
    1582729
    1583         diagout, diagoutd, diagoutvd = compute_clwvl(var0, qtotv, dnamesvar,dvnamesvar)
     730        diagout, diagoutd, diagoutvd = diag.compute_clwvl(var0, qtotv, dnamesvar,dvnamesvar)
    1584731
    1585732        # Removing the nonChecking variable-dimensions from the initial list
    … …  
    1639786        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
    1640787
    1641         diagout, diagoutd, diagoutvd = compute_prw(var0, var1, dnamesvar,dvnamesvar)
     788        diagout, diagoutd, diagoutvd = diag.compute_prw(var0, var1, dnamesvar,dvnamesvar)
    1642789
    1643790        # Removing the nonChecking variable-dimensions from the initial list
    … …  
    1667814        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
    1668815
    1669         diagout, diagoutd, diagoutvd = compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
     816        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
    1670817        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
    1671818
    … …  
    18771024# Global attributes
    18781025##
    1879 atvar = ncvar.set_attribute(newnc, 'program', 'diagnostics.py')
    1880 atvar = ncvar.set_attribute(newnc, 'version', '1.0')
    1881 atvar = ncvar.set_attribute(newnc, 'author', 'Fita Borrell, Lluis')
    1882 atvar = ncvar.set_attribute(newnc, 'institution', 'Laboratoire Meteorologie ' +      \
    1883   'Dynamique')
    1884 atvar = ncvar.set_attribute(newnc, 'university', 'Universite Pierre et Marie ' +     \
    1885   'Curie -- Jussieu')
    1886 atvar = ncvar.set_attribute(newnc, 'centre', 'Centre national de la recherche ' +    \
    1887   'scientifique')
    1888 atvar = ncvar.set_attribute(newnc, 'city', 'Paris')
    1889 atvar = ncvar.set_attribute(newnc, 'original_file', opts.ncfile)
     1026add_global_PyNCplot(newnc, main, None, '2.0')
    18901027
    18911028gorigattrs = ncobj.ncattrs()
    1892 
    18931029for attr in gorigattrs:
    18941030    attrv = ncobj.getncattr(attr)
Note: See TracChangeset for help on using the changeset viewer.