# Tools for the compute of diagnostics
# L. Fita, CIMA. CONICET-UBA, CNRS UMI-IFAECI, Buenos Aires, Argentina

# Available general pupose diagnostics (model independent) providing (varv1, varv2, ..., dimns, dimvns)
# compute_accum: Function to compute the accumulation of a variable
# compute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following 
#   newmicro.F90 from LMDZ compute_clt(cldfra, pres, dimns, dimvns)
# compute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ
# compute_clivi: Function to compute cloud-ice water path (clivi)
# compute_clwvl: Function to compute condensed water path (clwvl)
# compute_deaccum: Function to compute the deaccumulation of a variable
# compute_mslp: Function to compute mslp: mean sea level pressure following p_interp.F90 from WRF
# compute_OMEGAw: Function to transform OMEGA [Pas-1] to velocities [ms-1]
# compute_prw: Function to compute water vapour path (prw)
# compute_rh: Function to compute relative humidity following 'Tetens' equation (T,P) ...'
# compute_td: Function to compute the dew point temperature
# compute_turbulence: Function to compute the rubulence term of the Taylor's decomposition ...'
# C_diagnostic: Class to compute generic variables
# compute_wd: Function to compute the wind direction 3D
# compute_wds: Function to compute the wind direction
# compute_wss: Function to compute the wind speed
# compute_WRFhur: Function to compute WRF relative humidity following Teten's equation
# compute_WRFta: Function to compute WRF air temperature
# compute_WRFtd: Function to compute WRF dew-point air temperature
# compute_WRFua: Function to compute geographical rotated WRF x-wind
# compute_WRFva: Function to compute geographical rotated WRF y-wind
# compute_WRFuava: Function to compute geographical rotated WRF 3D winds
# compute_WRFuas: Function to compute geographical rotated WRF 2-meter x-wind
# compute_WRFvas: Function to compute geographical rotated WRF 2-meter y-wind
# compute_WRFuasvas: Fucntion to compute geographical rotated WRF 2-meter winds
# derivate_centered: Function to compute the centered derivate of a given field
# Forcompute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction 
#   following newmicro.F90 from LMDZ via Fortran subroutine
# Forcompute_clt: Function to compute the total cloud fraction following 
#   'newmicro.F90' from LMDZ via a Fortran module
# Forcompute_fog_K84: Computation of fog and visibility following Kunkel, (1984)
# Forcompute_fog_RUC: Computation of fog and visibility following RUC method Smirnova, (2000)
# Forcompute_fog_FRAML50: fog and visibility following Gultepe and Milbrandt, (2010)
# Forcompute_potevap_orPM: Function to compute potential evapotranspiration following
#   Penman-Monteith formulation implemented in ORCHIDEE
# Forcompute_psl_ptarget: Function to compute the sea-level pressure following 
#   target_pressure value found in `p_interp.F'
# Forcompute_zmla_gen: Function to compute the boundary layer height following a 
#   generic method with Fortran
# Forcompute_zwind: Function to compute the wind at a given height following the 
#   power law method
# Forcompute_zwind_log: Function to compute the wind at a given height following the 
#   logarithmic law method
# Forcompute_zwindMO: Function to compute the wind at a given height following the 
#   Monin-Obukhov theory
# W_diagnostic: Class to compute WRF diagnostics variables

# Others just providing variable values
# var_cllmh: Fcuntion to compute cllmh on a 1D column
# var_clt: Function to compute the total cloud fraction following 'newmicro.F90' from
#   LMDZ using 1D vertical column values
# var_convini: Function returns convective initialization (pr(t) > 0.0001) in time units
# var_hur: Function to compute relative humidity following 'August - Roche - Magnus' 
#   formula
# var_hur_Uhus: Function to compute relative humidity following 'August-Roche-Magnus'
#   formula using hus
# var_mslp: Fcuntion to compute mean sea-level pressure
# var_td: Function to compute dew-point air temperature from temperature and pressure
#   values
# var_td_Uhus: Function to compute dew-point air temperature from temperature and 
#   pressure values using hus
# var_timemax: This function returns the time at which variable reaches its maximum in time
#   units
# var_timeoverthres: This function returns the time at which (varv(t) > thres) in time units
# var_virtualTemp: This function returns virtual temperature in K, 
# var_WRFtime: Function to copmute CFtimes from WRFtime variable
# var_wd: Function to compute the wind direction
# var_wd: Function to compute the wind speed
# rotational_z: z-component of the rotatinoal of horizontal vectorial field
# turbulence_var: Function to compute the Taylor's decomposition turbulence term from
#   a a given variable

import numpy as np
from netCDF4 import Dataset as NetCDFFile
import os
import re
import nc_var_tools as ncvar
import generic_tools as gen
import datetime as dtime
import module_ForDiag as fdin
import module_ForDef as fdef

main = 'diag_tools.py'
errormsg = 'ERROR -- error -- ERROR -- error'
warnmsg = 'WARNING -- warning -- WARNING -- warning'

# Constants
grav = fdef.module_definitions.grav

# Available WRFiag
Wavailablediags = ['hur', 'p', 'ta', 'td', 'ua', 'va', 'uas', 'vas', 'wd', 'ws', 'zg']

# Available General diagnostics
Cavailablediags = ['hur', 'hur_Uhus', 'td', 'td_Uhus', 'wd', 'ws']

# Gneral information
##
def reduce_spaces(string):
    """ Function to give words of a line of text removing any extra space
    """
    values = string.replace('\n','').split(' ')
    vals = []
    for val in values:
         if len(val) > 0:
             vals.append(val)

    return vals

def variable_combo(varn,combofile):
    """ Function to provide variables combination from a given variable name
      varn= name of the variable
      combofile= ASCII file with the combination of variables
        [varn] [combo]
          [combo]: '@' separated list of variables to use to generate [varn]
            [WRFdt] to get WRF time-step (from general attributes)
    >>> variable_combo('WRFprls','/home/lluis/PY/diagnostics.inf')
    deaccum@RAINNC@XTIME@prnc
    """
    fname = 'variable_combo'

    if varn == 'h':
        print fname + '_____________________________________________________________'
        print variable_combo.__doc__
        quit()

    if not os.path.isfile(combofile):
        print errormsg
        print '  ' + fname + ": file with combinations '" + combofile +              \
          "' does not exist!!"
        quit(-1)

    objf = open(combofile, 'r')

    found = False
    for line in objf:
        linevals = reduce_spaces(line)
        varnf = linevals[0]
        combo = linevals[1].replace('\n','')
        if varn == varnf: 
            found = True
            break

    if not found:
        print errormsg
        print '  ' + fname + ": variable '" + varn + "' not found in '" + combofile +\
          "' !!"
        combo='ERROR'

    objf.close()

    return combo

# Mathematical operators
##
def compute_accum(varv, dimns, dimvns):
    """ Function to compute the accumulation of a variable
    compute_accum(varv, dimnames, dimvns)
      [varv]= values to accum (assuming [t,])
      [dimns]= list of the name of the dimensions of the [varv]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [varv]
    """
    fname = 'compute_accum'

    deacdims = dimns[:]
    deacvdims = dimvns[:]

    slicei = []
    slicee = []

    Ndims = len(varv.shape)
    for iid in range(0,Ndims):
        slicei.append(slice(0,varv.shape[iid]))
        slicee.append(slice(0,varv.shape[iid]))

    slicee[0] = np.arange(varv.shape[0])
    slicei[0] = np.arange(varv.shape[0])
    slicei[0][1:varv.shape[0]] = np.arange(varv.shape[0]-1)

    vari = varv[tuple(slicei)]
    vare = varv[tuple(slicee)]

    ac = vari*0.
    for it in range(1,varv.shape[0]):
        ac[it,] = ac[it-1,] + vare[it,]

    return ac, deacdims, deacvdims

def compute_deaccum(varv, dimns, dimvns):
    """ Function to compute the deaccumulation of a variable
    compute_deaccum(varv, dimnames, dimvns)
      [varv]= values to deaccum (assuming [t,])
      [dimns]= list of the name of the dimensions of the [varv]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [varv]
    """
    fname = 'compute_deaccum'

    deacdims = dimns[:]
    deacvdims = dimvns[:]

    slicei = []
    slicee = []

    Ndims = len(varv.shape)
    for iid in range(0,Ndims):
        slicei.append(slice(0,varv.shape[iid]))
        slicee.append(slice(0,varv.shape[iid]))

    slicee[0] = np.arange(varv.shape[0])
    slicei[0] = np.arange(varv.shape[0])
    slicei[0][1:varv.shape[0]] = np.arange(varv.shape[0]-1)

    vari = varv[tuple(slicei)]
    vare = varv[tuple(slicee)]

    deac = vare - vari

    return deac, deacdims, deacvdims

def derivate_centered(var,dim,dimv):
    """ Function to compute the centered derivate of a given field
      centered derivate(n) = (var(n-1) + var(n+1))/(2*dn).
    [var]= variable
    [dim]= which dimension to compute the derivate
    [dimv]= dimension values (can be of different dimension of [var])
    >>> derivate_centered(np.arange(16).reshape(4,4)*1.,1,1.)
    [[  0.   1.   2.   0.]
     [  0.   5.   6.   0.]
     [  0.   9.  10.   0.]
     [  0.  13.  14.   0.]]
    """

    fname = 'derivate_centered'
    
    vark = var.dtype

    if hasattr(dimv, "__len__"):
# Assuming that the last dimensions of var [..., N, M] are the same of dimv [N, M]
        if len(var.shape) != len(dimv.shape):
            dimvals = np.zeros((var.shape), dtype=vark)
            if len(var.shape) - len(dimv.shape) == 1:
                for iz in range(var.shape[0]):
                    dimvals[iz,] = dimv
            elif len(var.shape) - len(dimv.shape) == 2:
                for it in range(var.shape[0]):
                    for iz in range(var.shape[1]):
                        dimvals[it,iz,] = dimv
            else:
                print errormsg
                print '  ' + fname + ': dimension difference between variable',      \
                  var.shape,'and variable with dimension values',dimv.shape,         \
                  ' not ready !!!'
                quit(-1)
        else:
            dimvals = dimv
    else:
# dimension values are identical everywhere! 
# from: http://stackoverflow.com/questions/16807011/python-how-to-identify-if-a-variable-is-an-array-or-a-scalar    
        dimvals = np.ones((var.shape), dtype=vark)*dimv

    derivate = np.zeros((var.shape), dtype=vark)
    if dim > len(var.shape) - 1:
        print errormsg
        print '  ' + fname + ': dimension',dim,' too big for given variable of ' +   \
          'shape:', var.shape,'!!!'
        quit(-1)

    slicebef = []
    sliceaft = []
    sliceder = []

    for id in range(len(var.shape)):
        if id == dim:
            slicebef.append(slice(0,var.shape[id]-2))
            sliceaft.append(slice(2,var.shape[id]))
            sliceder.append(slice(1,var.shape[id]-1))
        else:
            slicebef.append(slice(0,var.shape[id]))
            sliceaft.append(slice(0,var.shape[id]))
            sliceder.append(slice(0,var.shape[id]))

    if hasattr(dimv, "__len__"):
        derivate[tuple(sliceder)] = (var[tuple(slicebef)] + var[tuple(sliceaft)])/   \
          ((dimvals[tuple(sliceaft)] - dimvals[tuple(slicebef)]))
        print (dimvals[tuple(sliceaft)] - dimvals[tuple(slicebef)])
    else:
        derivate[tuple(sliceder)] = (var[tuple(slicebef)] + var[tuple(sliceaft)])/   \
          (2.*dimv)

#    print 'before________'
#    print var[tuple(slicebef)]

#    print 'after________'
#    print var[tuple(sliceaft)]

    return derivate

def rotational_z(Vx,Vy,pos):
    """ z-component of the rotatinoal of horizontal vectorial field
    \/ x (Vx,Vy,Vz) = \/xVy - \/yVx
    [Vx]= Variable component x
    [Vy]=  Variable component y
    [pos]= poisition of the grid points
    >>> rotational_z(np.arange(16).reshape(4,4)*1., np.arange(16).reshape(4,4)*1., 1.)
    [[  0.   1.   2.   0.]
     [ -4.   0.   0.  -7.]
     [ -8.   0.   0. -11.]
     [  0.  13.  14.   0.]]
    """

    fname =  'rotational_z'

    ndims = len(Vx.shape)
    rot1 = derivate_centered(Vy,ndims-1,pos)
    rot2 = derivate_centered(Vx,ndims-2,pos)

    rot = rot1 - rot2

    return rot

# Diagnostics
##
def Forcompute_cape_afwa(ta, hur, pa, zg, hgt, parcelm, dimns, dimvns):
    """ Function to compute the CAPE, CIN, ZLFC, PLFC, LI following WRF 
      'phys/module_diaf_afwa.F' methodology
    Forcompute_cape_afwa(ta, hur, pa, hgt, zsfc, parcelm, dimns, dimvns)
      [ta]= air-temperature values (assuming [[t],z,y,x]) [K]
      [pa]= pressure values (assuming [[t],z,y,x]) [Pa]
      [zg]= gopotential height (assuming [[t],z,y,x]) [gpm]
      [hgt]= topographical height (assuming [m]
      [parcelm]= method of air-parcel to use
      [dimns]= list of the name of the dimensions of [pa]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [pa]
    """
    fname = 'Forcompute_cape_afwa'

    psldims = dimns[:]
    pslvdims = dimvns[:]

    if len(pa.shape) == 4:
        cape = np.zeros((pa.shape[0],pa.shape[2],pa.shape[3]), dtype=np.float)
        cin = np.zeros((pa.shape[0],pa.shape[2],pa.shape[3]), dtype=np.float)
        zlfc = np.zeros((pa.shape[0],pa.shape[2],pa.shape[3]), dtype=np.float)
        plfc = np.zeros((pa.shape[0],pa.shape[2],pa.shape[3]), dtype=np.float)
        li = np.zeros((pa.shape[0],pa.shape[2],pa.shape[3]), dtype=np.float)

        dx = pa.shape[3]
        dy = pa.shape[2]
        dz = pa.shape[1]
        dt = pa.shape[0]
        psldims.pop(1)
        pslvdims.pop(1)

        pcape,pcin,pzlfc,pplfc,pli= fdin.module_fordiagnostics.compute_cape_afwa4d(  \
          ta=ta[:].transpose(), hur=hur[:].transpose(), press=pa[:].transpose(),     \
          zg=zg[:].transpose(), hgt=hgt.transpose(), parcelmethod=parcelm, d1=dx,    \
          d2=dy, d3=dz, d4=dt)
        cape = pcape.transpose()
        cin = pcin.transpose()
        zlfc = pzlfc.transpose()
        plfc = pplfc.transpose()
        li = pli.transpose()
    else:
        print errormsg
        print '  ' + fname + ': rank', len(pa.shape), 'not ready !!'
        print '  it only computes 4D [t,z,y,x] rank values'
        quit(-1)

    return cape, cin, zlfc, plfc, li, psldims, pslvdims



def var_clt(cfra):
    """ Function to compute the total cloud fraction following 'newmicro.F90' from 
      LMDZ using 1D vertical column values
      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
    """
    ZEPSEC=1.0E-12

    fname = 'var_clt'

    zclear = 1.
    zcloud = 0.

    dz = cfra.shape[0]
    for iz in range(dz):
        zclear =zclear*(1.-np.max([cfra[iz],zcloud]))/(1.-np.min([zcloud,1.-ZEPSEC]))
        clt = 1. - zclear
        zcloud = cfra[iz]

    return clt

def compute_clt(cldfra, dimns, dimvns):
    """ Function to compute the total cloud fraction following 'newmicro.F90' from 
      LMDZ
    compute_clt(cldfra, dimnames)
      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
      [dimns]= list of the name of the dimensions of [cldfra]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [cldfra]
    """
    fname = 'compute_clt'

    cltdims = dimns[:]
    cltvdims = dimvns[:]

    if len(cldfra.shape) == 4:
        clt = np.zeros((cldfra.shape[0],cldfra.shape[2],cldfra.shape[3]),            \
          dtype=np.float)
        dx = cldfra.shape[3]
        dy = cldfra.shape[2]
        dz = cldfra.shape[1]
        dt = cldfra.shape[0]
        cltdims.pop(1)
        cltvdims.pop(1)

        for it in range(dt):
            for ix in range(dx):
                for iy in range(dy):
                    zclear = 1.
                    zcloud = 0.
                    gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
                    clt[it,iy,ix] = var_clt(cldfra[it,:,iy,ix])

    else:
        clt = np.zeros((cldfra.shape[1],cldfra.shape[2]), dtype=np.float)
        dx = cldfra.shape[2]
        dy = cldfra.shape[1]
        dy = cldfra.shape[0]
        cltdims.pop(0)
        cltvdims.pop(0)
        for ix in range(dx):
            for iy in range(dy):
                zclear = 1.
                zcloud = 0.
                gen.percendone(ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
                clt[iy,ix] = var_clt(cldfra[:,iy,ix])

    return clt, cltdims, cltvdims

def Forcompute_clt(cldfra, dimns, dimvns):
    """ Function to compute the total cloud fraction following 'newmicro.F90' from 
      LMDZ via a Fortran module
    compute_clt(cldfra, dimnames)
      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
      [dimns]= list of the name of the dimensions of [cldfra]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [cldfra]
    """
    fname = 'Forcompute_clt'

    cltdims = dimns[:]
    cltvdims = dimvns[:]


    if len(cldfra.shape) == 4:
        clt = np.zeros((cldfra.shape[0],cldfra.shape[2],cldfra.shape[3]),            \
          dtype=np.float)
        dx = cldfra.shape[3]
        dy = cldfra.shape[2]
        dz = cldfra.shape[1]
        dt = cldfra.shape[0]
        cltdims.pop(1)
        cltvdims.pop(1)

        clt = fdin.module_fordiagnostics.compute_clt4d2(cldfra[:])

    else:
        clt = np.zeros((cldfra.shape[1],cldfra.shape[2]), dtype=np.float)
        dx = cldfra.shape[2]
        dy = cldfra.shape[1]
        dy = cldfra.shape[0]
        cltdims.pop(0)
        cltvdims.pop(0)

        clt = fdin.module_fordiagnostics.compute_clt3d1(cldfra[:])

    return clt, cltdims, cltvdims

def var_cllmh(cfra, p):
    """ Fcuntion to compute cllmh on a 1D column
    """

    fname = 'var_cllmh'

    ZEPSEC =1.0E-12
    prmhc = 440.*100.
    prmlc = 680.*100.

    zclearl = 1.
    zcloudl = 0.
    zclearm = 1.
    zcloudm = 0.
    zclearh = 1.
    zcloudh = 0.

    dvz = cfra.shape[0]

    cllmh = np.ones((3), dtype=np.float)

    for iz in range(dvz):
        if p[iz] < prmhc:
            cllmh[2] = cllmh[2]*(1.-np.max([cfra[iz], zcloudh]))/(1.-                \
              np.min([zcloudh,1.-ZEPSEC]))
            zcloudh = cfra[iz]
        elif p[iz] >= prmhc and p[iz] < prmlc:
            cllmh[1] = cllmh[1]*(1.-np.max([cfra[iz], zcloudm]))/(1.-                \
              np.min([zcloudm,1.-ZEPSEC]))
            zcloudm = cfra[iz]
        elif p[iz] >= prmlc:
            cllmh[0] = cllmh[0]*(1.-np.max([cfra[iz], zcloudl]))/(1.-                \
              np.min([zcloudl,1.-ZEPSEC]))
            zcloudl = cfra[iz]

    cllmh = 1.- cllmh

    return cllmh

def Forcompute_cllmh(cldfra, pres, dimns, dimvns):
    """ Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ via Fortran subroutine
    compute_clt(cldfra, pres, dimns, dimvns)
      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
      [pres] = pressure field
      [dimns]= list of the name of the dimensions of [cldfra]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [cldfra]
    """
    fname = 'Forcompute_cllmh'

    cllmhdims = dimns[:]
    cllmhvdims = dimvns[:]

    if len(cldfra.shape) == 4:
        dx = cldfra.shape[3]
        dy = cldfra.shape[2]
        dz = cldfra.shape[1]
        dt = cldfra.shape[0]
        cllmhdims.pop(1)
        cllmhvdims.pop(1)

        cllmh = fdin.module_fordiagnostics.compute_cllmh4d2(cldfra[:], pres[:])
        
    else:
        dx = cldfra.shape[2]
        dy = cldfra.shape[1]
        dz = cldfra.shape[0]
        cllmhdims.pop(0)
        cllmhvdims.pop(0)

        cllmh = fdin.module_fordiagnostics.compute_cllmh3d1(cldfra[:], pres[:])

    return cllmh, cllmhdims, cllmhvdims

def compute_cllmh(cldfra, pres, dimns, dimvns):
    """ Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ
    compute_clt(cldfra, pres, dimns, dimvns)
      [cldfra]= cloud fraction values (assuming [[t],z,y,x])
      [pres] = pressure field
      [dimns]= list of the name of the dimensions of [cldfra]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [cldfra]
    """
    fname = 'compute_cllmh'

    cllmhdims = dimns[:]
    cllmhvdims = dimvns[:]

    if len(cldfra.shape) == 4:
        dx = cldfra.shape[3]
        dy = cldfra.shape[2]
        dz = cldfra.shape[1]
        dt = cldfra.shape[0]
        cllmhdims.pop(1)
        cllmhvdims.pop(1)

        cllmh = np.ones(tuple([3, dt, dy, dx]), dtype=np.float)

        for it in range(dt):
            for ix in range(dx):
                for iy in range(dy):
                    gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
                    cllmh[:,it,iy,ix] = var_cllmh(cldfra[it,:,iy,ix], pres[it,:,iy,ix])
        
    else:
        dx = cldfra.shape[2]
        dy = cldfra.shape[1]
        dz = cldfra.shape[0]
        cllmhdims.pop(0)
        cllmhvdims.pop(0)

        cllmh = np.ones(tuple([3, dy, dx]), dtype=np.float)

        for ix in range(dx):
            for iy in range(dy):
                gen.percendone(ix*dy + iy,dx*dy, 5, 'diagnosted')
                cllmh[:,iy,ix] = var_cllmh(cldfra[:,iy,ix], pres[:,iy,ix])

    return cllmh, cllmhdims, cllmhvdims

def compute_clivi(dens, qtot, dimns, dimvns):
    """ Function to compute cloud-ice water path (clivi)
      [dens] = density [in kgkg-1] (assuming [t],z,y,x)
      [qtot] = added mixing ratio of all cloud-ice species in [kgkg-1] (assuming [t],z,y,x)
      [dimns]= list of the name of the dimensions of [q]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [q]
    """
    fname = 'compute_clivi'

    clividims = dimns[:]
    clivivdims = dimvns[:]

    if len(qtot.shape) == 4:
        clividims.pop(1)
        clivivdims.pop(1)
    else:
        clividims.pop(0)
        clivivdims.pop(0)

    data1 = dens*qtot
    clivi = np.sum(data1, axis=1)

    return clivi, clividims, clivivdims


def compute_clwvl(dens, qtot, dimns, dimvns):
    """ Function to compute condensed water path (clwvl)
      [dens] = density [in kgkg-1] (assuming [t],z,y,x)
      [qtot] = added mixing ratio of all cloud-water species in [kgkg-1] (assuming [t],z,y,x)
      [dimns]= list of the name of the dimensions of [q]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [q]
    """
    fname = 'compute_clwvl'

    clwvldims = dimns[:]
    clwvlvdims = dimvns[:]

    if len(qtot.shape) == 4:
        clwvldims.pop(1)
        clwvlvdims.pop(1)
    else:
        clwvldims.pop(0)
        clwvlvdims.pop(0)

    data1 = dens*qtot
    clwvl = np.sum(data1, axis=1)

    return clwvl, clwvldims, clwvlvdims

def var_virtualTemp (temp,rmix):
    """ This function returns virtual temperature in K, 
      temp: temperature [K]
      rmix: mixing ratio in [kgkg-1]
    """

    fname = 'var_virtualTemp'

    virtual=temp*(0.622+rmix)/(0.622*(1.+rmix))

    return virtual

def var_convini(pr, time, dimns, dimvns):
    """ This function returns convective initialization (pr(t) > 0.0001) in time units
      pr: precipitation fux [kgm-2s-1]
      time: time in CF coordinates
    """
    fname = 'var_convini'

    dt = pr.shape[0]
    dy = pr.shape[1]
    dx = pr.shape[2]

    vardims = dimns[:]
    varvdims = dimvns[:]

    vardims.pop(0)
    varvdims.pop(0)

    prmin = 0.0001
    convini = np.ones((dy, dx), dtype=np.float)*gen.fillValueF
    for it in range(dt):
# NOT working ?
#        convini = np.where(convini == gen.fillValueF and pr[it,:,:] >= prmin,        \
#          time[it], fillValueF)
        for j in range(dy):
            for i in range(dx):
                if convini[j,i] == gen.fillValueF and pr[it,j,i] >= prmin:
                    convini[j,i] = time[it]
                    break

    return convini, vardims, varvdims

def var_timemax(varv, time, dimns, dimvns):
    """ This function returns the time at which variable reaches its maximum in time
        units
      varv: values of the variable to use
      time: time in CF coordinates
    """
    fname = 'var_timemax'

    dt = varv.shape[0]
    dy = varv.shape[1]
    dx = varv.shape[2]

    vardims = dimns[:]
    varvdims = dimvns[:]

    vardims.pop(0)
    varvdims.pop(0)

    timemax = np.ones((dy, dx), dtype=np.float)*gen.fillValueF
    varmax = np.max(varv, axis=0)
    for j in range(dy):
        for i in range(dx):
            it = gen.index_vec(varv[:,j,i], varmax[j,i])
            timemax[j,i] = time[it]

    return timemax, vardims, varvdims

def var_timeoverthres(varv, time, thres, dimns, dimvns):
    """ This function returns the time at which (varv(t) > thres) in time units
      varv: values of the variable to use
      time: time in CF coordinates
      thres: threshold to overpass
    """
    fname = 'var_timeoverthres'

    dt = varv.shape[0]
    dy = varv.shape[1]
    dx = varv.shape[2]

    vardims = dimns[:]
    varvdims = dimvns[:]

    vardims.pop(0)
    varvdims.pop(0)

    timeoverthres = np.ones((dy, dx), dtype=np.float)*gen.fillValueF
    for it in range(dt):
        for j in range(dy):
            for i in range(dx):
                if timeoverthres[j,i] == gen.fillValueF and varv[it,j,i] >= thres:
                    timeoverthres[j,i] = time[it]
                    break

    return timeoverthres, vardims, varvdims

def Forcompute_zint(var, zinterlev, zweights, dimns, dimvns):
    """ Function to compute a vertical integration of volumetric quantities
    Forcompute_mrso(smois, dsoil, dimns, dimvns)
      [var]= values (assuming [[t],z,y,x]) [volumetric units]
      [zinterlev]= depth of each layer (assuming [z]) [same z units as var]
      [zweights]= weights to apply to each level (just in case...)
      [dimns]= list of the name of the dimensions of [smois]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [smois]
    """
    fname = 'Forcompute_zint'

    vardims = dimns[:]
    varvdims = dimvns[:]

    if len(var.shape) == 4:
        zint = np.zeros((var.shape[0],var.shape[2],var.shape[3]), dtype=np.float)
        dx = var.shape[3]
        dy = var.shape[2]
        dz = var.shape[1]
        dt = var.shape[0]
        vardims.pop(1)
        varvdims.pop(1)

        zintvart=fdin.module_fordiagnostics.compute_zint4d(var4d=var[:].transpose(), \
          dlev=zinterlev[:].transpose(), zweight=zweights[:].transpose(), d1=dx,     \
          d2=dy, d3=dz, d4=dt)
        zintvar = zintvart.transpose()
    else:
        print errormsg
        print '  ' + fname + ': rank', len(var.shape), 'not ready !!'
        print '  it only computes 4D [t,z,y,x] rank values'
        quit(-1)

    return zintvar, vardims, varvdims

def var_mslp(pres, psfc, ter, tk, qv):
    """ Function to compute mslp on a 1D column
    """

    fname = 'var_mslp'

    N = 1.0
    expon=287.04*.0065/9.81
    pref = 40000.

# First find where about 400 hPa is located
    dz=len(pres) 

    kref = -1
    pinc = pres[0] - pres[dz-1]

    if pinc < 0.:
        for iz in range(1,dz):
            if pres[iz-1] >= pref and pres[iz] < pref: 
                kref = iz
                break
    else:
        for iz in range(dz-1):
            if pres[iz] >= pref and pres[iz+1] < pref: 
                kref = iz
                break

    if kref == -1:
        print errormsg
        print '  ' + fname + ': no reference pressure:',pref,'found!!'
        print '    values:',pres[:]
        quit(-1)

    mslp = 0.

# We are below both the ground and the lowest data level.

# First, find the model level that is closest to a "target" pressure
# level, where the "target" pressure is delta-p less that the local
# value of a horizontally smoothed surface pressure field.  We use
# delta-p = 150 hPa here. A standard lapse rate temperature profile
# passing through the temperature at this model level will be used
# to define the temperature profile below ground.  This is similar
# to the Benjamin and Miller (1990) method, using  
# 700 hPa everywhere for the "target" pressure.

# ptarget = psfc - 15000.
    ptarget = 70000.
    dpmin=1.e4
    kupper = 0
    if pinc > 0.:
        for iz in range(dz-1,0,-1):
            kupper = iz
            dp=np.abs( pres[iz] - ptarget )
            if dp < dpmin: exit
            dpmin = np.min([dpmin, dp])
    else:
        for iz in range(dz):
            kupper = iz
            dp=np.abs( pres[iz] - ptarget )
            if dp < dpmin: exit
            dpmin = np.min([dpmin, dp])

    pbot=np.max([pres[0], psfc])
#    zbot=0.

#    tbotextrap=tk(i,j,kupper,itt)*(pbot/pres_field(i,j,kupper,itt))**expon
#    tvbotextrap=virtual(tbotextrap,qv(i,j,1,itt))

#    data_out(i,j,itt,1) = (zbot+tvbotextrap/.0065*(1.-(interp_levels(1)/pbot)**expon))
    tbotextrap = tk[kupper]*(psfc/ptarget)**expon
    tvbotextrap = var_virtualTemp(tbotextrap, qv[kupper])
    mslp = psfc*( (tvbotextrap+0.0065*ter)/tvbotextrap)**(1./expon)

    return mslp

def compute_mslp(pressure, psurface, terrain, temperature, qvapor, dimns, dimvns):
    """ Function to compute mslp: mean sea level pressure following p_interp.F90 from WRF
    var_mslp(pres, ter, tk, qv, dimns, dimvns)
      [pressure]= pressure field [Pa] (assuming [[t],z,y,x])
      [psurface]= surface pressure field [Pa]
      [terrain]= topography [m]
      [temperature]= temperature [K]
      [qvapor]= water vapour mixing ratio [kgkg-1]
      [dimns]= list of the name of the dimensions of [cldfra]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [pres]
    """

    fname = 'compute_mslp'

    mslpdims = list(dimns[:])
    mslpvdims = list(dimvns[:])

    if len(pressure.shape) == 4:
        mslpdims.pop(1)
        mslpvdims.pop(1)
    else:
        mslpdims.pop(0)
        mslpvdims.pop(0)

    if len(pressure.shape) == 4:
        dx = pressure.shape[3]
        dy = pressure.shape[2]
        dz = pressure.shape[1]
        dt = pressure.shape[0]

        mslpv = np.zeros(tuple([dt, dy, dx]), dtype=np.float)

# Terrain... to 2D !
        terval = np.zeros(tuple([dy, dx]), dtype=np.float)
        if len(terrain.shape) == 3:
            terval = terrain[0,:,:]
        else:
            terval = terrain

        for ix in range(dx):
            for iy in range(dy):
                if terval[iy,ix] > 0.:
                    for it in range(dt):
                        mslpv[it,iy,ix] = var_mslp(pressure[it,:,iy,ix],             \
                          psurface[it,iy,ix], terval[iy,ix], temperature[it,:,iy,ix],\
                          qvapor[it,:,iy,ix])

                        gen.percendone(it*dx*dy + ix*dy + iy, dx*dy*dt, 5, 'diagnosted')
                else:
                    mslpv[:,iy,ix] = psurface[:,iy,ix]

    else:
        dx = pressure.shape[2]
        dy = pressure.shape[1]
        dz = pressure.shape[0]

        mslpv = np.zeros(tuple([dy, dx]), dtype=np.float)

# Terrain... to 2D !
        terval = np.zeros(tuple([dy, dx]), dtype=np.float)
        if len(terrain.shape) == 3:
            terval = terrain[0,:,:]
        else:
            terval = terrain

        for ix in range(dx):
            for iy in range(dy):
                gen.percendone(ix*dy + iy,dx*dy, 5, 'diagnosted')
                if terval[iy,ix] > 0.:
                    mslpv[iy,ix] = var_mslp(pressure[:,iy,ix], psurface[iy,ix],      \
                      terval[iy,ix], temperature[:,iy,ix], qvapor[:,iy,ix])
                else:
                    mslpv[iy,ix] = psfc[iy,ix]

    return mslpv, mslpdims, mslpvdims

def Forcompute_psl_ecmwf(ps, hgt, ta1, pa2, unpa1, dimns, dimvns):
    """ Function to compute the sea-level pressure following Mats Hamrud and Philippe Courtier [Pa]
    Forcompute_psl_ptarget(ps, hgt, ta1, pa2, unpa1, dimns, dimvns)
      [ps]= surface pressure values (assuming [[t],y,x]) [Pa]
      [hgt]= opography (assuming [y,x]) [m]
      [ta1]= air-temperature values at first half-mass level (assuming [[t],y,x]) [K]
      [pa2]= pressure values at second full-mass levels (assuming [[t],y,x]) [Pa]
      [unpa1]= pressure values at first half-mass levels (assuming [[t],y,x]) [Pa]
      [dimns]= list of the name of the dimensions of [pa]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [pa]
    """
    fname = 'Forcompute_psl_ecmwf'

    vardims = dimns[:]
    varvdims = dimvns[:]

    if len(pa2.shape) == 3:
        psl = np.zeros((pa2.shape[0],pa2.shape[1],pa2.shape[2]), dtype=np.float)
        dx = pa2.shape[2]
        dy = pa2.shape[1]
        dt = pa2.shape[0]
        pslt= fdin.module_fordiagnostics.compute_psl_ecmwf( ps=ps[:].transpose(),    \
          hgt=hgt[:].transpose(), t=ta1[:].transpose(), press=pa2[:].transpose(),    \
          unpress=unpa1[:].transpose(), d1=dx, d2=dy, d4=dt)
        psl = pslt.transpose()
    else:
        print errormsg
        print '  ' + fname + ': rank', len(pa2.shape), 'not ready !!'
        print '  it only computes 3D [t,y,x] rank values'
        quit(-1)

    return psl, vardims, varvdims

def Forcompute_psl_ptarget(pa, ps, ta, hgt, qv, target_pressure, dimns, dimvns):
    """ Function to compute the sea-level pressure following target_pressure value
      found in `p_interp.F'
    Forcompute_psl_ptarget(pa, ps, ta, hgt, qv, dimns, dimvns)
      [pa]= pressure values (assuming [[t],z,y,x]) [Pa]
      [ps]= surface pressure values (assuming [[t],y,x]) [Pa]
      [ta]= air-temperature values (assuming [[t],z,y,x]) [K]
      [hgt]= opography (assuming [y,x]) [m]
      [qv]= water vapour mixing ratio (assuming [[t],z,y,x]) [kgkg-1]
      [dimns]= list of the name of the dimensions of [pa]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [pa]
    """
    fname = 'Forcompute_psl_ptarget'

    psldims = dimns[:]
    pslvdims = dimvns[:]

    if len(pa.shape) == 4:
        psl = np.zeros((pa.shape[0],pa.shape[2],pa.shape[3]), dtype=np.float)
        dx = pa.shape[3]
        dy = pa.shape[2]
        dz = pa.shape[1]
        dt = pa.shape[0]
        psldims.pop(1)
        pslvdims.pop(1)

        pslt= fdin.module_fordiagnostics.compute_psl_ptarget4d2(                     \
          press=pa[:].transpose(), ps=ps[:].transpose(), hgt=hgt[:].transpose(),     \
          ta=ta[:].transpose(), qv=qv[:].transpose(), ptarget=target_pressure,       \
          d1=dx, d2=dy, d3=dz, d4=dt)
        psl = pslt.transpose()
    else:
        print errormsg
        print '  ' + fname + ': rank', len(pa.shape), 'not ready !!'
        print '  it only computes 4D [t,z,y,x] rank values'
        quit(-1)

    return psl, psldims, pslvdims

def Forcompute_zmla_gen(theta, qratio, zpl, hgt, dimns, dimvns):
    """ Function to compute the boundary layer height following a generic method 
      with Fortran
    Forcompute_zmla_gen(theta, qratio, zpl, hgt, zmla, dimns, dimvns)
      [theta]= potential air-temperature values (assuming [[t],z,y,x]) [K]
      [qratio]= water mixing ratio (assuming [[t],z,y,x]) [kgkg-1]
      [zpl]= height from sea level (assuming [[t],z,y,x]) [m]
      [hgt]= topographical height (assuming [m]
      [zmla]= boundary layer height [m]
      [dimns]= list of the name of the dimensions of [theta]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [theta]
    """
    fname = 'Forcompute_zmla_gen'

    zmladims = dimns[:]
    zmlavdims = dimvns[:]

    if len(theta.shape) == 4:
        zmla= np.zeros((theta.shape[0],theta.shape[2],theta.shape[3]), dtype=np.float)

        dx = theta.shape[3]
        dy = theta.shape[2]
        dz = theta.shape[1]
        dt = theta.shape[0]
        zmladims.pop(1)
        zmlavdims.pop(1)

        pzmla= fdin.module_fordiagnostics.compute_zmla_generic4d(                    \
          tpot=theta[:].transpose(), qratio=qratio[:].transpose(),                   \
          z=zpl[:].transpose(), hgt=hgt.transpose(), d1=dx, d2=dy, d3=dz, d4=dt)
        zmla = pzmla.transpose()
    else:
        print errormsg
        print '  ' + fname + ': rank', len(theta.shape), 'not ready !!'
        print '  it only computes 4D [t,z,y,x] rank values'
        quit(-1)

    return zmla, zmladims, zmlavdims

def Forcompute_zwind(ua, va, z, uas, vas, sina, cosa, zval, dimns, dimvns):
    """ Function to compute the wind at a given height following the power law method
    Forcompute_zwind(ua, va, zsl, uas, vas, hgt, sina, cosa, zval, dimns, dimvns)
      [ua]= x-component of unstaggered 3D wind (assuming [[t],z,y,x]) [ms-1]
      [va]= y-component of unstaggered 3D wind (assuming [[t],z,y,x]) [ms-1]
      [z]= height above surface [m]
      [uas]= x-component of unstaggered 10 m wind (assuming [[t],z,y,x]) [ms-1]
      [vas]= y-component of unstaggered 10 m wind (assuming [[t],z,y,x]) [ms-1]
      [sina]= local sine of map rotation [1.]
      [cosa]= local cosine of map rotation [1.]
      [zval]= desired height for winds [m]
      [dimns]= list of the name of the dimensions of [ua]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [ua]
    """
    fname = 'Forcompute_zwind'

    vardims = dimns[:]
    varvdims = dimvns[:]

    if len(ua.shape) == 4:
        var1= np.zeros((ua.shape[0],ua.shape[2],ua.shape[3]), dtype=np.float)
        var2= np.zeros((ua.shape[0],ua.shape[2],ua.shape[3]), dtype=np.float)

        dx = ua.shape[3]
        dy = ua.shape[2]
        dz = ua.shape[1]
        dt = ua.shape[0]
        vardims.pop(1)
        varvdims.pop(1)

        pvar1, pvar2= fdin.module_fordiagnostics.compute_zwind4d(ua=ua.transpose(),  \
          va=va[:].transpose(), z=z[:].transpose(), uas=uas.transpose(),             \
          vas=vas.transpose(), sina=sina.transpose(), cosa=cosa.transpose(),         \
          zextrap=zval, d1=dx, d2=dy, d3=dz, d4=dt)
        var1 = pvar1.transpose()
        var2 = pvar2.transpose()
    else:
        print errormsg
        print '  ' + fname + ': rank', len(ua.shape), 'not ready !!'
        print '  it only computes 4D [t,z,y,x] rank values'
        quit(-1)

    return var1, var2, vardims, varvdims

def Forcompute_zwind_log(ua, va, z, uas, vas, sina, cosa, zval, dimns, dimvns):
    """ Function to compute the wind at a given height following the logarithmic law method
    Forcompute_zwind(ua, va, zsl, uas, vas, hgt, sina, cosa, zval, dimns, dimvns)
      [ua]= x-component of unstaggered 3D wind (assuming [[t],z,y,x]) [ms-1]
      [va]= y-component of unstaggered 3D wind (assuming [[t],z,y,x]) [ms-1]
      [z]= height above surface [m]
      [uas]= x-component of unstaggered 10 m wind (assuming [[t],z,y,x]) [ms-1]
      [vas]= y-component of unstaggered 10 m wind (assuming [[t],z,y,x]) [ms-1]
      [sina]= local sine of map rotation [1.]
      [cosa]= local cosine of map rotation [1.]
      [zval]= desired height for winds [m]
      [dimns]= list of the name of the dimensions of [ua]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [ua]
    """
    fname = 'Forcompute_zwind_log'

    vardims = dimns[:]
    varvdims = dimvns[:]

    if len(ua.shape) == 4:
        var1= np.zeros((ua.shape[0],ua.shape[2],ua.shape[3]), dtype=np.float)
        var2= np.zeros((ua.shape[0],ua.shape[2],ua.shape[3]), dtype=np.float)

        dx = ua.shape[3]
        dy = ua.shape[2]
        dz = ua.shape[1]
        dt = ua.shape[0]
        vardims.pop(1)
        varvdims.pop(1)

        pvar1, pvar2= fdin.module_fordiagnostics.compute_zwind_log4d(                \
          ua=ua.transpose(), va=va[:].transpose(), z=z[:].transpose(),               \
          uas=uas.transpose(), vas=vas.transpose(), sina=sina.transpose(),           \
          cosa=cosa.transpose(), zextrap=zval, d1=dx, d2=dy, d3=dz, d4=dt)
        var1 = pvar1.transpose()
        var2 = pvar2.transpose()
    else:
        print errormsg
        print '  ' + fname + ': rank', len(ua.shape), 'not ready !!'
        print '  it only computes 4D [t,z,y,x] rank values'
        quit(-1)

    return var1, var2, vardims, varvdims

def Forcompute_zwindMO(ust, znt, rmol, uas, vas, sina, cosa, zval, dimns, dimvns):
    """ Function to compute the wind at a given height following the Monin-Obukhov theory
    Forcompute_zwind(ust, znt, rmol, uas, vas, sina, cosa, zval, dimns, dimvns)
      [ust]: u* in similarity theory (assuming [[t],y,x]) [ms-1]
      [znt]: thermal time-varying roughness length (assuming [[t],y,x]) [m]
      [rmol]: inverse of Obukhov length (assuming [[t],y,x]) [m-1]
      [uas]= x-component of unstaggered 10 m wind (assuming [[t],y,x]) [ms-1]
      [vas]= y-component of unstaggered 10 m wind (assuming [[t],y,x]) [ms-1]
      [sina]= local sine of map rotation [1.]
      [cosa]= local cosine of map rotation [1.]
      [zval]= desired height for winds [m]
      [dimns]= list of the name of the dimensions of [uas]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [uas]
    """
    fname = 'Forcompute_zwindMO'

    vardims = dimns[:]
    varvdims = dimvns[:]

    if len(uas.shape) == 3:
        var1= np.zeros((uas.shape[0],uas.shape[1],uas.shape[2]), dtype=np.float)
        var2= np.zeros((uas.shape[0],uas.shape[1],uas.shape[2]), dtype=np.float)

        dx = uas.shape[2]
        dy = uas.shape[1]
        dt = uas.shape[0]

        pvar1, pvar2 = fdin.module_fordiagnostics.compute_zwindmo3d(                 \
          ust=ust.transpose(), znt=znt[:].transpose(), rmol=rmol[:].transpose(),     \
          uas=uas.transpose(), vas=vas.transpose(), sina=sina.transpose(),           \
          cosa=cosa.transpose(), newz=zval, d1=dx, d2=dy, d3=dt)
        var1 = pvar1.transpose()
        var2 = pvar2.transpose()
    else:
        print errormsg
        print '  ' + fname + ': rank', len(uas.shape), 'not ready !!'
        print '  it only computes 3D [t,y,x] rank values'
        quit(-1)

    return var1, var2, vardims, varvdims

def Forcompute_potevap_orPM(rho1, ust, uas, vas, tas, ps, qv1, dimns, dimvns):
    """ Function to compute potential evapotranspiration following Penman-Monteith 
      formulation implemented in ORCHIDEE in src_sechiba/enerbil.f90
    Forcompute_potevap_orPM(rho1, uas, vas, tas, ps, qv2, qv1, dimns, dimvns)
      [rho1]= air-density at the first layer (assuming [[t],y,m]) [kgm-3]
      [ust]= u* in similarity theory (assuming [[t],y,x]) [ms-1]
      [uas]= x-component of unstaggered 10 m wind (assuming [[t],y,x]) [ms-1]
      [vas]= y-component of unstaggered 10 m wind (assuming [[t],y,x]) [ms-1]
      [tas]= 2m air temperature [K]
      [ps]= surface pressure [Pa]
      [qv1]= mixing ratio at the first atmospheric layer [kgkg-1]
      [dimns]= list of the name of the dimensions of [uas]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [uas]
    """
    fname = 'Forcompute_potevap_orPM'

    vardims = dimns[:]
    varvdims = dimvns[:]

    if len(uas.shape) == 3:
        var1= np.zeros((uas.shape[0],uas.shape[1],uas.shape[2]), dtype=np.float)
        var2= np.zeros((uas.shape[0],uas.shape[1],uas.shape[2]), dtype=np.float)

        dx = uas.shape[2]
        dy = uas.shape[1]
        dt = uas.shape[0]

        pvar = fdin.module_fordiagnostics.compute_potevap_orpm3d(                    \
          rho1=rho1.transpose(), ust=ust.transpose(), uas=uas.transpose(),           \
          vas=vas.transpose(), tas=tas.transpose(), ps=ps.transpose(),               \
          qv1=qv1.transpose(), d1=dx, d2=dy, d3=dt)
        var = pvar.transpose()
    else:
        print errormsg
        print '  ' + fname + ': rank', len(uas.shape), 'not ready !!'
        print '  it only computes 3D [t,y,x] rank values'
        quit(-1)

    return var, vardims, varvdims

def Forcompute_fog_K84(qcloud, qice, dimns, dimvns):
    """ Function to compute fog and visibility following Kunkel, (1984)
    Forcompute_fog_K84(qcloud, qice, dimns, dimvns)
      [qcloud]= cloud mixing ratio [kgk-1]
      [qice]= ice mixing ratio [kgk-1]
      [dimns]= list of the name of the dimensions of [uas]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [qcloud]
    """
    fname = 'Forcompute_fog_K84'

    vardims = dimns[:]
    varvdims = dimvns[:]

    if len(qcloud.shape) == 4:
        var= np.zeros((qcloud.shape[0],qcloud.shape[2],qcloud.shape[3]), dtype=np.float)

        dx = qcloud.shape[3]
        dy = qcloud.shape[2]
        dz = qcloud.shape[1]
        dt = qcloud.shape[0]
        vardims.pop(1)
        varvdims.pop(1)

        pvar1, pvar2 = fdin.module_fordiagnostics.compute_fog_k84(                   \
          qc=qcloud[:,0,:,:].transpose(), qi=qice[:,0,:,:].transpose(), d1=dx, d2=dy,\
          d3=dt)
        var1 = pvar1.transpose()
        var2 = pvar2.transpose()
    else:
        print errormsg
        print '  ' + fname + ': rank', len(qcloud.shape), 'not ready !!'
        print '  it only computes 4D [t,z,y,x] rank values'
        quit(-1)

    return var1, var2, vardims, varvdims

def Forcompute_fog_RUC(qvapor, temp, pres, dimns, dimvns):
    """ Function to compute fog and visibility following RUC method Smirnova, (2000)
    Forcompute_fog_RUC(qcloud, qice, dimns, dimvns)
      [qvapor]= water vapor mixing ratio [kgk-1]
      [temp]= temperature [K]
      [pres]= pressure [Pa]
      [dimns]= list of the name of the dimensions of [uas]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [qcloud]
    """
    fname = 'Forcompute_fog_RUC'

    vardims = dimns[:]
    varvdims = dimvns[:]

    if len(qvapor.shape) == 4:
        var= np.zeros((qvapor.shape[0],qvapor.shape[2],qvapor.shape[3]), dtype=np.float)

        dx = qvapor.shape[3]
        dy = qvapor.shape[2]
        dz = qvapor.shape[1]
        dt = qvapor.shape[0]
        vardims.pop(1)
        varvdims.pop(1)

        pvar1, pvar2 = fdin.module_fordiagnostics.compute_fog_ruc(                   \
          qv=qvapor[:,0,:,:].transpose(), ta=temp[:,0,:,:].transpose(),              \
          pres=pres[:,0,:,:].transpose(), d1=dx, d2=dy, d3=dt)
        var1 = pvar1.transpose()
        var2 = pvar2.transpose()
    elif len(qvapor.shape) == 3:
        var= np.zeros((qvapor.shape[0],qvapor.shape[1],qvapor.shape[2]), dtype=np.float)

        dx = qvapor.shape[2]
        dy = qvapor.shape[1]
        dt = qvapor.shape[0]

        pvar1, pvar2 = fdin.module_fordiagnostics.compute_fog_ruc(                   \
          qv=qvapor[:].transpose(), ta=temp[:].transpose(), pres=pres[:].transpose(),\
          d1=dx, d2=dy, d3=dt)
        var1 = pvar1.transpose()
        var2 = pvar2.transpose()
    else:
        print errormsg
        print '  ' + fname + ': rank', len(qcloud.shape), 'not ready !!'
        print '  it only computes 4D [t,z,y,x] or 3D [t,z,y,x] rank values'
        quit(-1)

    return var1, var2, vardims, varvdims

def Forcompute_fog_FRAML50(qvapor, temp, pres, dimns, dimvns):
    """ Function to compute fog (vis < 1km) and visibility following FRAM-L 50 % prob
         Gultepe, and Milbrandt, (2010), J. Appl. Meteor. Climatol.
    Forcompute_fog_FRAML50(qvapor, temp, pres, dimns, dimvns)
      [qvapor]= vapor mixing ratio [kgk-1]
      [temp]= temperature [K]
      [pres]= pressure [Pa]
      [dimns]= list of the name of the dimensions of [uas]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [qvapor]
    """
    fname = 'Forcompute_fog_FRAML50'

    vardims = dimns[:]
    varvdims = dimvns[:]

    if len(qvapor.shape) == 4:
        var= np.zeros((qvapor.shape[0],qvapor.shape[2],qvapor.shape[3]), dtype=np.float)

        dx = qvapor.shape[3]
        dy = qvapor.shape[2]
        dz = qvapor.shape[1]
        dt = qvapor.shape[0]
        vardims.pop(1)
        varvdims.pop(1)

        pvar1, pvar2 = fdin.module_fordiagnostics.compute_fog_framl50(               \
          qv=qvapor[:,0,:,:].transpose(), ta=temp[:,0,:,:].transpose(),              \
          pres=pres[:,0,:,:].transpose(), d1=dx, d2=dy, d3=dt)
        var1 = pvar1.transpose()
        var2 = pvar2.transpose()
    elif len(qvapor.shape) == 3:
        var= np.zeros((qvapor.shape[0],qvapor.shape[1],qvapor.shape[2]), dtype=np.float)

        dx = qvapor.shape[2]
        dy = qvapor.shape[1]
        dt = qvapor.shape[0]

        pvar1, pvar2 = fdin.module_fordiagnostics.compute_fog_framl50(               \
          qv=qvapor[:].transpose(), ta=temp[:].transpose(), pres=pres[:].transpose(),\
          d1=dx, d2=dy, d3=dt)
        var1 = pvar1.transpose()
        var2 = pvar2.transpose()
    else:
        print errormsg
        print '  ' + fname + ': rank', len(qvapor.shape), 'not ready !!'
        print '  it only computes 4D [t,z,y,x] or 3D [t,y,x] rank values'
        quit(-1)

    return var1, var2, vardims, varvdims

####### ###### ##### #### ### ## # END Fortran diagnostics

def compute_OMEGAw(omega, p, t, dimns, dimvns):
    """ Function to transform OMEGA [Pas-1] to velocities [ms-1]
      tacking: https://www.ncl.ucar.edu/Document/Functions/Contributed/omega_to_w.shtml
      [omega] = vertical velocity [in ms-1] (assuming [t],z,y,x)
      [p] = pressure in [Pa] (assuming [t],z,y,x)
      [t] = temperature in [K] (assuming [t],z,y,x)
      [dimns]= list of the name of the dimensions of [q]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [q]
    """
    fname = 'compute_OMEGAw'

    rgas = 287.058            # J/(kg-K) => m2/(s2 K)
    g    = 9.80665            # m/s2

    wdims = dimns[:]
    wvdims = dimvns[:]

    rho  = p/(rgas*t)         # density => kg/m3
    w    = -omega/(rho*g)     

    return w, wdims, wvdims

def compute_prw(dens, q, dimns, dimvns):
    """ Function to compute water vapour path (prw)
      [dens] = density [in kgkg-1] (assuming [t],z,y,x)
      [q] = mixing ratio in [kgkg-1] (assuming [t],z,y,x)
      [dimns]= list of the name of the dimensions of [q]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [q]
    """
    fname = 'compute_prw'

    prwdims = dimns[:]
    prwvdims = dimvns[:]

    if len(q.shape) == 4:
        prwdims.pop(1)
        prwvdims.pop(1)
    else:
        prwdims.pop(0)
        prwvdims.pop(0)

    data1 = dens*q
    prw = np.sum(data1, axis=1)

    return prw, prwdims, prwvdims

def compute_rh(p, t, q, dimns, dimvns):
    """ Function to compute relative humidity following 'Tetens' equation (T,P) ...'
      [t]= temperature (assuming [[t],z,y,x] in [K])
      [p] = pressure field (assuming in [hPa])
      [q] = mixing ratio in [kgkg-1]
      [dimns]= list of the name of the dimensions of [t]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [t]
    """
    fname = 'compute_rh'

    rhdims = dimns[:]
    rhvdims = dimvns[:]

    data1 = 10.*0.6112*np.exp(17.67*(t-273.16)/(t-29.65))
    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)

    rh = q/data2

    return rh, rhdims, rhvdims

def compute_td(p, temp, qv, dimns, dimvns):
    """ Function to compute the dew point temperature
      [p]= pressure [Pa]
      [temp]= temperature [C]
      [qv]= mixing ratio [kgkg-1]
      [dimns]= list of the name of the dimensions of [p]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [p]
    """
    fname = 'compute_td'

#    print '    ' + fname + ': computing dew-point temperature from TS as t and Tetens...'
# tacking from: http://en.wikipedia.org/wiki/Dew_point
    tk = temp
    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)

    rh = qv/data2
                
    pa = rh * data1
    td = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))

    tddims = dimns[:]
    tdvdims = dimvns[:]

    return td, tddims, tdvdims

def var_WRFtime(timewrfv, refdate='19491201000000', tunitsval='minutes'):
    """ Function to copmute CFtimes from WRFtime variable
      refdate= [YYYYMMDDMIHHSS] format of reference date
      tunitsval= CF time units
      timewrfv= matrix string values of WRF 'Times' variable
    """
    fname = 'var_WRFtime'

    yrref=refdate[0:4]
    monref=refdate[4:6]
    dayref=refdate[6:8]
    horref=refdate[8:10]
    minref=refdate[10:12]
    secref=refdate[12:14]

    refdateS = yrref + '-' + monref + '-' + dayref + ' ' + horref + ':' + minref +   \
      ':' + secref

    dt = timewrfv.shape[0]
    WRFtime = np.zeros((dt), dtype=np.float)

    for it in range(dt):
        wrfdates = gen.datetimeStr_conversion(timewrfv[it,:],'WRFdatetime', 'matYmdHMS')
        WRFtime[it] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)

    tunits = tunitsval + ' since ' + refdateS

    return WRFtime, tunits

def turbulence_var(varv, dimvn, dimn):
    """ Function to compute the Taylor's decomposition turbulence term from a a given variable
      x*=<x^2>_t-(<X>_t)^2
    turbulence_var(varv,dimn)
      varv= values of the variable
      dimvn= names of the dimension of the variable
      dimn= names of the dimensions (as a dictionary with 'X', 'Y', 'Z', 'T')
    >>> turbulence_var(np.arange((27)).reshape(3,3,3),['time','y','x'],{'T':'time', 'Y':'y', 'X':'x'})
    [[ 54.  54.  54.]
     [ 54.  54.  54.]
     [ 54.  54.  54.]]
    """
    fname = 'turbulence_varv'

    timedimid = dimvn.index(dimn['T'])

    varv2 = varv*varv

    vartmean = np.mean(varv, axis=timedimid)
    var2tmean = np.mean(varv2, axis=timedimid)

    varvturb = var2tmean - (vartmean*vartmean)

    return varvturb

def compute_turbulence(v, dimns, dimvns):
    """ Function to compute the rubulence term of the Taylor's decomposition ...'
      x*=<x^2>_t-(<X>_t)^2
      [v]= variable (assuming [[t],z,y,x])
      [dimns]= list of the name of the dimensions of [v]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [v]
    """
    fname = 'compute_turbulence'

    turbdims = dimns[:]
    turbvdims = dimvns[:]

    turbdims.pop(0)
    turbvdims.pop(0)

    v2 = v*v

    vartmean = np.mean(v, axis=0)
    var2tmean = np.mean(v2, axis=0)

    turb = var2tmean - (vartmean*vartmean)

    return turb, turbdims, turbvdims

def compute_wd(u, v, dimns, dimvns):
    """ Function to compute the wind direction
      [u]= W-E wind direction [ms-1, knot, ...]
      [v]= N-S wind direction [ms-1, knot, ...]
      [dimns]= list of the name of the dimensions of [u]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [u]
    """
    fname = 'compute_wds'

#    print '    ' + fname + ': computing wind direction as ATAN2(v,u) ...'
    theta = np.arctan2(v,u)
    theta = np.where(theta < 0., theta + 2.*np.pi, theta)

    var = 360.*theta/(2.*np.pi)

    vardims = dimns[:]
    varvdims = dimvns[:]

    return var, vardims, varvdims

def compute_wds(u, v, dimns, dimvns):
    """ Function to compute the wind direction
      [u]= W-E wind direction [ms-1, knot, ...]
      [v]= N-S wind direction [ms-1, knot, ...]
      [dimns]= list of the name of the dimensions of [u]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [u]
    """
    fname = 'compute_wds'

#    print '    ' + fname + ': computing wind direction as ATAN2(v,u) ...'
    theta = np.arctan2(v,u)
    theta = np.where(theta < 0., theta + 2.*np.pi, theta)

    wds = 360.*theta/(2.*np.pi)

    wdsdims = dimns[:]
    wdsvdims = dimvns[:]

    return wds, wdsdims, wdsvdims

def compute_wss(u, v, dimns, dimvns):
    """ Function to compute the wind speed
      [u]= W-E wind direction [ms-1, knot, ...]
      [v]= N-S wind direction [ms-1, knot, ...]
      [dimns]= list of the name of the dimensions of [u]
      [dimvns]= list of the name of the variables with the values of the 
        dimensions of [u]
    """
    fname = 'compute_wss'

#    print '    ' + fname + ': computing wind speed as SQRT(v**2 + u**2) ...'
    wss = np.sqrt(u*u + v*v)

    wssdims = dimns[:]
    wssvdims = dimvns[:]

    return wss, wssdims, wssvdims

def timeunits_seconds(dtu):
    """ Function to transform a time units to seconds
    timeunits_seconds(timeuv)
      [dtu]= time units value to transform in seconds
    """
    fname='timunits_seconds'

    if dtu == 'years':
        times = 365.*24.*3600.
    elif dtu == 'weeks':
        times = 7.*24.*3600.
    elif dtu == 'days':
        times = 24.*3600.
    elif dtu == 'hours':
        times = 3600.
    elif dtu == 'minutes':
        times = 60.
    elif dtu == 'seconds':
        times = 1.
    elif dtu == 'miliseconds':
        times = 1./1000.
    else:
        print errormsg
        print '  ' + fname  + ": time units '" + dtu + "' not ready !!"
        quit(-1)

    return times

def compute_WRFhur(t, p, qv, dimns, dimvns):
    """ Function to compute WRF relative humidity following Teten's equation
      t= orginal WRF temperature
      p= original WRF pressure (P + PB)
      formula:
          temp = theta*(p/p0)**(R/Cp)

    """
    fname = 'compute_WRFtd'

    tk = (t+300.)*(p/fdef.module_definitions.p0ref)**(fdef.module_definitions.rcp)
    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)

    rh = qv/data2
                
    dnamesvar = ['Time','bottom_top','south_north','west_east']
    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)

    return rh, dnamesvar, dvnamesvar

def compute_WRFua(u, v, sina, cosa, dimns, dimvns):
    """ Function to compute geographical rotated WRF 3D winds
      u= orginal WRF x-wind
      v= orginal WRF y-wind
      sina= original WRF local sinus of map rotation
      cosa= original WRF local cosinus of map rotation
      formula:
        ua = u*cosa-va*sina
        va = u*sina+va*cosa
    """
    fname = 'compute_WRFua'

    var0 = u
    var1 = v
    var2 = sina
    var3 = cosa

    # un-staggering variables
    unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
    ua = np.zeros(tuple(unstgdims), dtype=np.float)
    unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
    unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
    unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
    unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])

    for iz in range(var0.shape[1]):
        ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2

    dnamesvar = ['Time','bottom_top','south_north','west_east']
    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)

    return ua, dnamesvar, dvnamesvar

def compute_WRFva(u, v, sina, cosa, dimns, dimvns):
    """ Function to compute geographical rotated WRF 3D winds
      u= orginal WRF x-wind
      v= orginal WRF y-wind
      sina= original WRF local sinus of map rotation
      cosa= original WRF local cosinus of map rotation
      formula:
        ua = u*cosa-va*sina
        va = u*sina+va*cosa
    """
    fname = 'compute_WRFva'

    var0 = u
    var1 = v
    var2 = sina
    var3 = cosa

    # un-staggering variables
    unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
    va = np.zeros(tuple(unstgdims), dtype=np.float)
    unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
    unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
    unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
    unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])

    for iz in range(var0.shape[1]):
        va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3

    dnamesvar = ['Time','bottom_top','south_north','west_east']
    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)

    return va, dnamesvar, dvnamesvar

def compute_WRFuava(u, v, sina, cosa, dimns, dimvns):
    """ Function to compute geographical rotated WRF 3D winds
      u= orginal WRF x-wind
      v= orginal WRF y-wind
      sina= original WRF local sinus of map rotation
      cosa= original WRF local cosinus of map rotation
      formula:
        ua = u*cosa-va*sina
        va = u*sina+va*cosa
    """
    fname = 'compute_WRFuava'

    var0 = u
    var1 = v
    var2 = sina
    var3 = cosa

    # un-staggering variables
    unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
    ua = np.zeros(tuple(unstgdims), dtype=np.float)
    va = np.zeros(tuple(unstgdims), dtype=np.float)
    unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
    unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
    unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
    unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])

    for iz in range(var0.shape[1]):
        ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
        va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3

    dnamesvar = ['Time','bottom_top','south_north','west_east']
    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)

    return ua, va, dnamesvar, dvnamesvar

def compute_WRFuas(u10, v10, sina, cosa, dimns, dimvns):
    """ Function to compute geographical rotated WRF 2-meter x-wind
      u10= orginal WRF 10m x-wind
      v10= orginal WRF 10m y-wind
      sina= original WRF local sinus of map rotation
      cosa= original WRF local cosinus of map rotation
      formula:
        uas = u10*cosa-va10*sina
        vas = u10*sina+va10*cosa
    """
    fname = 'compute_WRFuas'

    var0 = u10
    var1 = v10
    var2 = sina
    var3 = cosa

    uas = np.zeros(var0.shape, dtype=np.float)
    vas = np.zeros(var0.shape, dtype=np.float)

    uas = var0*var3 - var1*var2

    dnamesvar = ['Time','south_north','west_east']
    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)

    return uas, dnamesvar, dvnamesvar

def compute_WRFvas(u10, v10, sina, cosa, dimns, dimvns):
    """ Function to compute geographical rotated WRF 2-meter y-wind
      u10= orginal WRF 10m x-wind
      v10= orginal WRF 10m y-wind
      sina= original WRF local sinus of map rotation
      cosa= original WRF local cosinus of map rotation
      formula:
        uas = u10*cosa-va10*sina
        vas = u10*sina+va10*cosa
    """
    fname = 'compute_WRFvas'

    var0 = u10
    var1 = v10
    var2 = sina
    var3 = cosa

    uas = np.zeros(var0.shape, dtype=np.float)
    vas = np.zeros(var0.shape, dtype=np.float)

    vas = var0*var2 + var1*var3

    dnamesvar = ['Time','south_north','west_east']
    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)

    return vas, dnamesvar, dvnamesvar

def compute_WRFuasvas(u10, v10, sina, cosa, dimns, dimvns):
    """ Function to compute geographical rotated WRF 2-meter winds
      u10= orginal WRF 10m x-wind
      v10= orginal WRF 10m y-wind
      sina= original WRF local sinus of map rotation
      cosa= original WRF local cosinus of map rotation
      formula:
        uas = u10*cosa-va10*sina
        vas = u10*sina+va10*cosa
    """
    fname = 'compute_WRFuasvas'

    var0 = u10
    var1 = v10
    var2 = sina
    var3 = cosa

    uas = np.zeros(var0.shape, dtype=np.float)
    vas = np.zeros(var0.shape, dtype=np.float)

    uas = var0*var3 - var1*var2
    vas = var0*var2 + var1*var3

    dnamesvar = ['Time','south_north','west_east']
    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)

    return uas, vas, dnamesvar, dvnamesvar

def compute_WRFta(t, p, dimns, dimvns):
    """ Function to compute WRF air temperature
      t= orginal WRF temperature
      p= original WRF pressure (P + PB)
      formula:
          temp = theta*(p/p0)**(R/Cp)

    """
    fname = 'compute_WRFta'

    ta = (t+300.)*(p/fdef.module_definitions.p0ref)**(fdef.module_definitions.rcp)

    dnamesvar = ['Time','bottom_top','south_north','west_east']
    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)

    return ta, dnamesvar, dvnamesvar

def compute_WRFtd(t, p, qv, dimns, dimvns):
    """ Function to compute WRF dew-point air temperature
      t= orginal WRF temperature
      p= original WRF pressure (P + PB)
      formula:
          temp = theta*(p/p0)**(R/Cp)

    """
    fname = 'compute_WRFtd'

    tk = (t+300.)*(p/fdef.module_definitions.p0ref)**(fdef.module_definitions.rcp)
    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)

    rh = qv/data2
                
    pa = rh * data1
    td = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))

    dnamesvar = ['Time','bottom_top','south_north','west_east']
    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)

    return td, dnamesvar, dvnamesvar

def compute_WRFwd(u, v, sina, cosa, dimns, dimvns):
    """ Function to compute the wind direction
      u= W-E wind direction [ms-1]
      v= N-S wind direction [ms-1]
      sina= original WRF local sinus of map rotation
      cosa= original WRF local cosinus of map rotation
    """
    fname = 'compute_WRFwd'
    var0 = u
    var1 = v
    var2 = sina
    var3 = cosa

    # un-staggering variables
    unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
    ua = np.zeros(tuple(unstgdims), dtype=np.float)
    va = np.zeros(tuple(unstgdims), dtype=np.float)
    unstgvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
    unstgvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
    unstgvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
    unstgvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])

    for iz in range(var0.shape[1]):
        ua[:,iz,:,:] = unstgvar0[:,iz,:,:]*var3 - unstgvar1[:,iz,:,:]*var2
        va[:,iz,:,:] = unstgvar0[:,iz,:,:]*var2 + unstgvar1[:,iz,:,:]*var3

    theta = np.arctan2(va,ua)
    theta = np.where(theta < 0., theta + 2.*np.pi, theta)

    wd = 360.*theta/(2.*np.pi)

    dnamesvar = ['Time','bottom_top','south_north','west_east']
    dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dimns,dimvns)

    return wd

####### Variables (as they arrive without dimensions staff)

def var_hur(p, t, q):
    """ Function to compute relative humidity following 'August - Roche - Magnus' formula
      [t]= temperature (assuming [[t],z,y,x] in [K])
      [p] = pressure field (assuming in [Pa])
      [q] = mixing ratio in [kgkg-1]
      ref.: M. G. Lawrence, BAMS, 2005, 225
    >>> print var_rh(101300., 290., 3.)
    0.250250256174
    """
    fname = 'var_hur'

    # Enthalpy of vaporization [Jkg-1]
    L = 2.453*10.**6
    # Gas constant for water vapor [JK-1Kg-1]
    Rv = 461.5 

    # Computing saturated pressure 
    #es = 10.*0.6112*np.exp(17.67*(t-273.16)/(t-29.65))
    #es = 6.11*np.exp(L/Rv*(1.-273./t)/273.)

    # August - Roche - Magnus formula
    es = 6.1094*np.exp(17.625*(t-273.15)/((t-273.15)+243.04))

    # Saturated mixing ratio [g/kg]
    ws = 0.622*es/(0.01*p-es)

    hur = q/(ws*1000.)

    #print 'q:', q[5,5,5,5], 't:', t[5,5,5,5], 'p:', p[5,5,5,5]
    #print 'es:', es[5,5,5,5], 'ws:', ws[5,5,5,5], 'hur:', hur[5,5,5,5]

    return hur

def var_hur_Uhus(p, t, hus):
    """ Function to compute relative humidity following 'August - Roche - Magnus' formula using hus
      [t]= temperature (assuming [[t],z,y,x] in [K])
      [p] = pressure field (assuming in [Pa])
      [hus] = specific humidty [1]
      ref.: M. G. Lawrence, BAMS, 2005, 225
    >>> print var_rh(101300., 290., 3.)
    0.250250256174
    """
    fname = 'var_hur_Uhus'

    # Enthalpy of vaporization [Jkg-1]
    L = 2.453*10.**6
    # Gas constant for water vapor [JK-1Kg-1]
    Rv = 461.5 

    # Computing saturated pressure 
    #es = 10.*0.6112*np.exp(17.67*(t-273.16)/(t-29.65))
    #es = 6.11*np.exp(L/Rv*(1.-273./t)/273.)

    # August - Roche - Magnus formula
    es = 6.1094*np.exp(17.625*(t-273.15)/((t-273.15)+243.04))

    # Saturated mixing ratio [g/kg]
    ws = 0.622*es/(0.01*p-es)

    # Mixing ratio
    q = hus/(1.+hus)

    hur = q/ws

    #print 'q:', q[5,5,5,5], 't:', t[5,5,5,5], 'p:', p[5,5,5,5]
    #print 'es:', es[5,5,5,5], 'ws:', ws[5,5,5,5], 'hur:', hur[5,5,5,5]

    return hur

def var_td(t, p, qv):
    """ Function to compute dew-point air temperature from temperature and pressure values
      t= temperature [K]
      p= pressure (Pa)
      formula:
          temp = theta*(p/p0)**(R/Cp)

    """
    fname = 'compute_td'

    tk = (t)*(p/fdef.module_definitions.p0ref)**(fdef.module_definitions.rcp)
    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)

    rh = qv/data2
                
    pa = rh * data1
    td = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))

    return td

def var_td_Uhus(t, p, hus):
    """ Function to compute dew-point air temperature from temperature and pressure values using hus
      t= temperature [K]
      hus= specific humidity [1]
      formula:
          temp = theta*(p/p0)**(R/Cp)

    """
    fname = 'compute_td'

    tk = (t)*(p/fdef.module_definitions.p0ref)**(fdef.module_definitions.rcp)
    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
    data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)

    qv = hus/(1.+hus)

    rh = qv/data2
                
    pa = rh * data1
    td = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))

    return td

def var_wd(u, v):
    """ Function to compute the wind direction
      [u]= W-E wind direction [ms-1, knot, ...]
      [v]= N-S wind direction [ms-1, knot, ...]
    """
    fname = 'var_wd'

    theta = np.arctan2(v,u)
    theta = np.where(theta < 0., theta + 2.*np.pi, theta)

    wd = 360.*theta/(2.*np.pi)

    return wd

def var_ws(u, v):
    """ Function to compute the wind speed
      [u]= W-E wind direction [ms-1, knot, ...]
      [v]= N-S wind direction [ms-1, knot, ...]
    """
    fname = 'var_ws'

    ws = np.sqrt(u*u + v*v)

    return ws

class C_diagnostic(object):
    """ Class to compute generic variables
      Cdiag: name of the diagnostic to compute
      ncobj: netcdf object with data
      sfcvars: dictionary with CF equivalencies of surface variables inside file
      vars3D: dictionary with CF equivalencies of 3D variables inside file
      dictdims: dictionary with CF equivalencies of dimensions inside file
        self.values = Values of the diagnostic
        self.dims = Dimensions of the diagnostic
        self.units = units of the diagnostic
        self.incvars = list of variables from the input netCDF object
    """ 
    def __init__(self, Cdiag, ncobj, sfcvars, vars3D, dictdims):
        fname = 'C_diagnostic'
        self.values = None
        self.dims = None
        self.incvars = ncobj.variables
        self.units = None

        if Cdiag == 'hur':
            """ Computing relative humidity
            """
            vn = 'hur'
            CF3Dvars = ['ta', 'plev', 'q']
            for v3D in CF3Dvars:
                if not vars3D.has_key(v3D):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + v3D +              \
                      "' attribution to compute '" + vn + "' !!"
                    print '  Equivalence of 3D variables provided _______'
                    gen.printing_dictionary(vars3D)
                    quit(-1)
                if not self.incvars.has_key(vars3D[v3D]):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + vars3D[v3D] +      \
                      "' in input file to compute '" + vn + "' !!"
                    print '  available variables:', self.incvars.keys()
                    print '  looking for variables _______'  
                    gen.printing_dictionary(vars3D)
                    quit(-1)

            ta = ncobj.variables[vars3D['ta']][:]
            p = ncobj.variables[vars3D['plev']][:]
            q = ncobj.variables[vars3D['q']][:]

            if len(ta.shape) != len(p.shape):
                p = gen.fill_Narray(p, ta*0., filldim=[0,2,3])

            self.values = var_hur(p, ta, q)
            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
              dictdims['lon']]
            self.units = '1'

        if Cdiag == 'hur_Uhus':
            """ Computing relative humidity using hus
            """
            vn = 'hur'
            CF3Dvars = ['ta', 'plev', 'hus']
            for v3D in CF3Dvars:
                if not vars3D.has_key(v3D):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + v3D +              \
                      "' attribution to compute '" + vn + "' !!"
                    print '  Equivalence of 3D variables provided _______'
                    gen.printing_dictionary(vars3D)
                    quit(-1)
                if not self.incvars.has_key(vars3D[v3D]):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + vars3D[v3D] +      \
                      "' in input file to compute '" + vn + "' !!"
                    print '  available variables:', self.incvars.keys()
                    print '  looking for variables _______'  
                    gen.printing_dictionary(vars3D)
                    quit(-1)

            ta = ncobj.variables[vars3D['ta']][:]
            p = ncobj.variables[vars3D['plev']][:]
            hus = ncobj.variables[vars3D['hus']][:]

            if len(ta.shape) != len(p.shape):
                p = gen.fill_Narray(p, ta*0., filldim=[0,2,3])

            self.values = var_hur_Uhus(p, ta, hus)
            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
              dictdims['lon']]
            self.units = '1'

        elif Cdiag == 'td':
            """ Computing dew-point temperature 
            """
            vn = 'td'
            CF3Dvars = ['ta', 'plev', 'q']
            for v3D in CF3Dvars:
                if not vars3D.has_key(v3D):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + v3D +              \
                      "' attribution to compute '" + vn + "' !!"
                    print '  Equivalence of 3D variables provided _______'
                    gen.printing_dictionary(vars3D)
                    quit(-1)
                if not self.incvars.has_key(vars3D[v3D]):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + vars3D[v3D] +      \
                      "' in input file to compute '" + vn + "' !!"
                    print '  available variables:', self.incvars.keys()
                    print '  looking for variables _______'  
                    gen.printing_dictionary(vars3D)
                    quit(-1)

            ta = ncobj.variables[vars3D['ta']][:]
            p = ncobj.variables[vars3D['plev']][:]
            q = ncobj.variables[vars3D['q']][:]

            if len(ta.shape) != len(p.shape):
                p = gen.fill_Narray(p, ta*0., filldim=[0,2,3])

            self.values = var_td(ta, p, q)
            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
              dictdims['lon']]
            self.units = 'K'

        elif Cdiag == 'td_Uhus':
            """ Computing dew-point temperature 
            """
            vn = 'td'
            CF3Dvars = ['ta', 'plev', 'hus']
            for v3D in CF3Dvars:
                if not vars3D.has_key(v3D):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + v3D +              \
                      "' attribution to compute '" + vn + "' !!"
                    print '  Equivalence of 3D variables provided _______'
                    gen.printing_dictionary(vars3D)
                    quit(-1)
                if not self.incvars.has_key(vars3D[v3D]):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + vars3D[v3D] +      \
                      "' in input file to compute '" + vn + "' !!"
                    print '  available variables:', self.incvars.keys()
                    print '  looking for variables _______'  
                    gen.printing_dictionary(vars3D)
                    quit(-1)

            ta = ncobj.variables[vars3D['ta']][:]
            p = ncobj.variables[vars3D['plev']][:]
            hus = ncobj.variables[vars3D['hus']][:]

            if len(ta.shape) != len(p.shape):
                p = gen.fill_Narray(p, ta*0., filldim=[0,2,3])

            self.values = var_td_Uhus(ta, p, hus)
            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
              dictdims['lon']]
            self.units = 'K'

        elif Cdiag == 'wd':
            """ Computing wind direction
            """
            vn = 'wd'
            CF3Dvars = ['ua', 'va']
            for v3D in CF3Dvars:
                if not vars3D.has_key(v3D):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + v3D +              \
                      "self.' attribution to compute '" + vn + "' !!"
                    print '  Equivalence of 3D variables provided _______'
                    gen.printing_dictionary(vars3D)
                    quit(-1)
                if not self.incvars.has_key(vars3D[v3D]):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + vars3D[v3D] +      \
                      "' in input file to compute '" + vn + "' !!"
                    print '  available variables:', self.incvars.keys()
                    print '  looking for variables _______'  
                    gen.printing_dictionary(vars3D)
                    quit(-1)
    
            ua = ncobj.variables[vars3D['ua']][:]
            va = ncobj.variables[vars3D['va']][:]
    
            self.values = var_wd(ua, va)
            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
              dictdims['lon']]
            self.units = 'degree'

        elif Cdiag == 'ws':
            """ Computing wind speed
            """
            vn = 'ws'
            CF3Dvars = ['ua', 'va']
            for v3D in CF3Dvars:
                if not vars3D.has_key(v3D):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + v3D +              \
                      "' attribution to compute '" + vn + "' !!"
                    print '  Equivalence of 3D variables provided _______'
                    gen.printing_dictionary(vars3D)
                    quit(-1)
                if not self.incvars.has_key(vars3D[v3D]):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + vars3D[v3D] +      \
                      "' in input file to compute '" + vn + "' !!"
                    print '  available variables:', self.incvars.keys()
                    print '  looking for variables _______'  
                    gen.printing_dictionary(vars3D)
                    quit(-1)

            ua = ncobj.variables[vars3D['ua']][:]
            va = ncobj.variables[vars3D['va']][:]

            self.values = var_ws(ua, va)
            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
              dictdims['lon']]
            self.units = ncobj.variables[vars3D['ua']].units

        else:
            print gen.errormsg
            print '  ' + fname + ": variable '" + Wdiag + "' not ready !!"
            print '  available ones:', Cavailablediags
            quit(-1)

class W_diagnostic(object):
    """ Class to compute WRF diagnostics variables
      Wdiag: name of the diagnostic to compute
      ncobj: netcdf object with data
      sfcvars: dictionary with CF equivalencies of surface variables inside file
      vars3D: dictionary with CF equivalencies of 3D variables inside file
      indims: list of dimensions inside file
      invardims: list of dimension-variables inside file
      dictdims: dictionary with CF equivalencies of dimensions inside file
        self.values = Values of the diagnostic
        self.dims = Dimensions of the diagnostic
        self.units = units of the diagnostic
        self.incvars = list of variables from the input netCDF object
    """    
    def __init__(self, Wdiag, ncobj, sfcvars, vars3D, indims, invardims, dictdims):
        fname = 'W_diagnostic'

        self.values = None
        self.dims = None
        self.incvars = ncobj.variables
        self.units = None

        if Wdiag == 'hur':
            """ Computing relative humidity
            """
            vn = 'hur'
            CF3Dvars = ['ta', 'hus']
            for v3D in CF3Dvars:
                if not vars3D.has_key(v3D):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + v3D +              \
                      "' attribution to compute '" + vn + "' !!"
                    print '  Equivalence of 3D variables provided _______'
                    gen.printing_dictionary(vars3D)
                    quit(-1)

            ta = ncobj.variables['T'][:]
            p = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
            hur = ncobj.variables['QVAPOR'][:]
    
            vals, dims, vdims = compute_WRFhur(ta, p, hur, indims, invardims)
            self.values = vals
            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
              dictdims['lon']]
            self.units = '1'

        elif Wdiag == 'p':
            """ Computing air pressure
            """
            vn = 'p'

            self.values = ncobj.variables['PB'][:] + ncobj.variables['P'][:]
            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
              dictdims['lon']]
            self.units = ncobj.variables['PB'].units

        elif Wdiag == 'ta':
            """ Computing air temperature 
            """
            vn = 'ta'
            CF3Dvars = ['ta']
            for v3D in CF3Dvars:
                if not vars3D.has_key(v3D):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + v3D +              \
                      "' attribution to compute '" + vn + "' !!"
                    print '  Equivalence of 3D variables provided _______'
                    gen.printing_dictionary(vars3D)
                    quit(-1)

            ta = ncobj.variables['T'][:]
            p = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
    
            vals, dims, vdims = compute_WRFta(ta, p, indims, invardims)
            self.values = vals
            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
              dictdims['lon']]
            self.units = 'K'

        elif Wdiag == 'td':
            """ Computing dew-point temperature 
            """
            vn = 'td'
            CF3Dvars = ['ta', 'hus']
            for v3D in CF3Dvars:
                if not vars3D.has_key(v3D):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + v3D +              \
                      "' attribution to compute '" + vn + "' !!"
                    print '  Equivalence of 3D variables provided _______'
                    gen.printing_dictionary(vars3D)
                    quit(-1)

            ta = ncobj.variables['T'][:]
            p = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
            hus = ncobj.variables['QVAPOR'][:]
    
            vals, dims, vdims = compute_WRFtd(ta, p, hus, indims, invardims)
            self.values = vals
            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
              dictdims['lon']]
            self.units = 'K'
    
        elif Wdiag == 'ua':
            """ Computing x-wind
            """
            vn = 'ua'
            CF3Dvars = ['ua', 'va']
            for v3D in CF3Dvars:
                if not vars3D.has_key(v3D):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + v3D +              \
                      "' attribution to compute '" + vn + "' !!"
                    print '  Equivalence of 3D variables provided _______'
                    gen.printing_dictionary(vars3D)
                    quit(-1)

            ua = ncobj.variables['U'][:]
            va = ncobj.variables['V'][:]
            sina = ncobj.variables['SINALPHA'][:]
            cosa = ncobj.variables['COSALPHA'][:]
    
            vals, dims, vdims = compute_WRFua(ua, va, sina, cosa, indims, invardims)
            self.values = vals
            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
              dictdims['lon']]
            self.units = ncobj.variables['U'].units
    
        elif Wdiag == 'uas':
            """ Computing 10m x-wind
            """
            vn = 'uas'
            CFsfcvars = ['uas', 'vas']
            for vsf in CFsfcvars:
                if not sfcvars.has_key(vsf):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + vsf +              \
                      "' attribution to compute '" + vn + "' !!"
                    print '  Equivalence of sfc variables provided _______'
                    gen.printing_dictionary(sfcvars)
                    quit(-1)
    
            uas = ncobj.variables['U10'][:]
            vas = ncobj.variables['V10'][:]
            sina = ncobj.variables['SINALPHA'][:]
            cosa = ncobj.variables['COSALPHA'][:]
    
            vals,dims,vdims = compute_WRFuas(uas, vas, sina, cosa, indims, invardims)
            self.values = vals
            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
              dictdims['lon']]
            self.units = ncobj.variables['U10'].units
    
        elif Wdiag == 'va':
            """ Computing y-wind
            """
            vn = 'ua'
            CF3Dvars = ['ua', 'va']
            for v3D in CF3Dvars:
                if not vars3D.has_key(v3D):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + v3D +              \
                      "' attribution to compute '" + vn + "' !!"
                    print '  Equivalence of 3D variables provided _______'
                    gen.printing_dictionary(vars3D)
                    quit(-1)
    
            ua = ncobj.variables['U'][:]
            va = ncobj.variables['V'][:]
            sina = ncobj.variables['SINALPHA'][:]
            cosa = ncobj.variables['COSALPHA'][:]
    
            vals, dims, vdims = compute_WRFva(ua, va, sina, cosa, indims, invardims)
            self.values = vals
            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
              dictdims['lon']]
            self.units = ncobj.variables['U'].units
    
        elif Wdiag == 'vas':
            """ Computing 10m y-wind
            """
            vn = 'uas'
            CFsfcvars = ['uas', 'vas']
            for vsf in CFsfcvars:
                if not sfcvars.has_key(vsf):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + vsf +              \
                      "' attribution to compute '" + vn + "' !!"
                    print '  Equivalence of sfc variables provided _______'
                    gen.printing_dictionary(sfcvars)
                    quit(-1)
    
            uas = ncobj.variables['U10'][:]
            vas = ncobj.variables['V10'][:]
            sina = ncobj.variables['SINALPHA'][:]
            cosa = ncobj.variables['COSALPHA'][:]
    
            vals,dims,vdims = compute_WRFvas(uas, vas, sina, cosa, indims, invardims)
            self.values = vals
            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
              dictdims['lon']]
            self.units = ncobj.variables['U10'].units

        elif Wdiag == 'wd':
            """ Computing wind direction
            """
            vn = 'wd'
            CF3Dvars = ['ua', 'va']
            for v3D in CF3Dvars:
                if not vars3D.has_key(v3D):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + v3D +              \
                      "' attribution to compute '" + vn + "' !!"
                    print '  Equivalence of 3D variables provided _______'
                    gen.printing_dictionary(vars3D)
                    quit(-1)

            ua = ncobj.variables['U10'][:]
            va = ncobj.variables['V10'][:]
            sina = ncobj.variables['SINALPHA'][:]
            cosa = ncobj.variables['COSALPHA'][:]
    
            vals, dims, vdims = compute_WRFwd(ua, va, sina, cosa, indims, invardims)
            self.values = vals
            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
              dictdims['lon']]
            self.units = 'degree'

        elif Wdiag == 'ws':
            """ Computing wind speed
            """
            vn = 'ws'
            CF3Dvars = ['ua', 'va']
            for v3D in CF3Dvars:
                if not vars3D.has_key(v3D):
                    print gen.errormsg
                    print '  ' + fname + ": missing variable '" + v3D +              \
                      "' attribution to compute '" + vn + "' !!"
                    print '  Equivalence of 3D variables provided _______'
                    gen.printing_dictionary(vars3D)
                    quit(-1)
    
            ua = ncobj.variables['U10'][:]
            va = ncobj.variables['V10'][:]
    
            self.values = var_ws(ua, va)
            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
              dictdims['lon']]
            self.units = ncobj.variables['U10'].units

        elif Wdiag == 'zg':
            """ Computing geopotential
            """
            vn = 'zg'

            self.values = ncobj.variables['PHB'][:] + ncobj.variables['PH'][:]
            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
              dictdims['lon']]
            self.units = ncobj.variables['PHB'].units

        else:
            print gen.errormsg
            print '  ' + fname + ": variable '" + Wdiag + "' not ready !!"
            print '  available ones:', Wavailablediags
            quit(-1)