# Python script to comput diagnostics
# L. Fita, LMD. CNR, UPMC-Jussieu, Paris, France
# File diagnostics.inf provides the combination of variables to get the desired diagnostic
#   To be used with module_ForDiagnostics.F90, module_ForDiagnosticsVars.F90, module_generic.F90
#      foudre: f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90 >& run_f2py.log
#      ciclad: f2py --f90flags="-fPIC" --f90exec=/usr/bin/gfortran -L/opt/canopy-1.3.0/Canopy_64bit/System/lib/ -L/usr/lib64/ -L/opt/canopy-1.3.0/Canopy_64bit/System/lib/ -m module_ForDiagnostics -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90 >& run_f2py.log

## e.g. # diagnostics.py -d 'Time@WRFtime,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
## e.g. # diagnostics.py -f /home/lluis/PY/diagnostics.inf -d variable_combo -v WRFprc

# 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 ...'
# compute_wds: Function to compute the wind direction
# compute_wss: Function to compute the wind speed
# compute_WRFuava: Function to compute geographical rotated WRF 3D winds
# compute_WRFuasvas: Fucntion to compute geographical rotated WRF 2-meter winds
# derivate_centered: Function to compute the centered derivate of a given field
# def 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_psl_ptarget: Function to compute the sea-level pressure following target_pressure value found in `p_interp.F'

# 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_mslp: Fcuntion to compute mean sea-level pressure
# var_virtualTemp: This function returns virtual temperature in K, 
# var_WRFtime: Function to copmute CFtimes from WRFtime variable
# 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

from optparse import OptionParser
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 diag_tools as diag

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

# Constants
grav = 9.81


####### ###### ##### #### ### ## #
comboinf="\nIF -d 'variable_combo', provides information of the combination to obtain -v [varn] with the ASCII file with the combinations as -f [combofile]"

parser = OptionParser()
parser.add_option("-f", "--netCDF_file", dest="ncfile", help="file to use", metavar="FILE")
parser.add_option("-d", "--dimensions", dest="dimns",  
  help="[dimtn]@[dtvn],[dimzn]@[dzvn],[...,[dimxn]@[dxvn]], ',' list with the couples [dimDn]@[dDvn], [dimDn], name of the dimension D and name of the variable [dDvn] with the values of the dimension ('WRFtime', for WRF time copmutation). NOTE: same order as in file!!!!" + comboinf, 
  metavar="LABELS")
parser.add_option("-v", "--variables", dest="varns", 
  help=" [varn1]|[var11]@[...[varN1]],[...,[varnM]|[var1M]@[...[varLM]]] ',' list of variables to compute [varnK] and its necessary ones [var1K]...[varPK]", metavar="VALUES")

(opts, args) = parser.parse_args()

#######    #######
## MAIN
    #######
availdiags = ['ACRAINTOT', 'accum', 'clt', 'cllmh', 'deaccum', 'LMDZrh', 'mslp',     \
  'OMEGAw', 'RAINTOT',                                                               \
  'rvors', 'td', 'turbulence', 'WRFcape_afwa', 'WRFclivi', 'WRFclwvl', 'WRFgeop',    \
  'WRFmrso', 'WRFp', 'WRFpsl_ecmwf',                                                 \
  'WRFpsl_ptarget', 'WRFrvors', 'WRFslw', 'ws', 'wds', 'wss', 'WRFheight',           \
  'WRFheightrel', 'WRFua', 'WRFva', 'WRFzwind', 'WRFzwind_log', 'WRFzwindMO']

methods = ['accum', 'deaccum']

# Variables not to check
NONcheckingvars = ['cllmh', 'deaccum', 'TSrhs', 'TStd', 'TSwds', 'TSwss', 'WRFbils', \
  'WRFclivi', 'WRFclwvl', 'WRFdens', 'WRFgeop',                                      \
  'WRFp', 'WRFtd',                                                                   \
  'WRFpos', 'WRFprc', 'WRFprls', 'WRFrh', 'LMDZrh', 'LMDZrhs', 'WRFrhs', 'WRFrvors', \
  'WRFt', 'WRFtime', 'WRFua', 'WRFva', 'WRFwds', 'WRFwss', 'WRFheight', 'WRFz']

NONchkvardims = ['WRFtime']

ofile = 'diagnostics.nc'

dimns = opts.dimns
varns = opts.varns

# Special method. knowing variable combination
##
if opts.dimns == 'variable_combo':
    print warnmsg
    print '  ' + main + ': knowing variable combination !!!'
    combination = variable_combo(opts.varns,opts.ncfile)
    print '     COMBO: ' + combination
    quit(-1)

if not os.path.isfile(opts.ncfile):
    print errormsg
    print '   ' + main + ": file '" + opts.ncfile + "' does not exist !!"
    quit(-1)

ncobj = NetCDFFile(opts.ncfile, 'r')

# Looking for specific variables that might be use in more than one diagnostic
WRFgeop_compute = False
WRFp_compute = False
WRFt_compute = False
WRFrh_compute = False
WRFght_compute = False
WRFdens_compute = False
WRFpos_compute = False
WRFtime_compute = False
WRFz_compute = False

# File creation
newnc = NetCDFFile(ofile,'w')

# dimensions
dimvalues = dimns.split(',')
dnames = []
dvnames = []

for dimval in dimvalues:
    dn = dimval.split('@')[0]
    dnv = dimval.split('@')[1]
    dnames.append(dn)
    dvnames.append(dnv)
    # Is there any dimension-variable which should be computed?
    if dnv == 'WRFgeop':WRFgeop_compute = True
    if dnv == 'WRFp': WRFp_compute = True
    if dnv == 'WRFt': WRFt_compute = True
    if dnv == 'WRFrh': WRFrh_compute = True
    if dnv == 'WRFght': WRFght_compute = True
    if dnv == 'WRFdens': WRFdens_compute = True
    if dnv == 'WRFpos': WRFpos_compute = True
    if dnv == 'WRFtime': WRFtime_compute = True
    if dnv == 'WRFz':WRFz_compute = True

# diagnostics to compute
diags = varns.split(',')
Ndiags = len(diags)

for idiag in range(Ndiags):
    if diags[idiag].split('|')[1].find('@') == -1:
        depvars = diags[idiag].split('|')[1]
        if depvars == 'WRFgeop':WRFgeop_compute = True
        if depvars == 'WRFp': WRFp_compute = True
        if depvars == 'WRFt': WRFt_compute = True
        if depvars == 'WRFrh': WRFrh_compute = True
        if depvars == 'WRFght': WRFght_compute = True
        if depvars == 'WRFdens': WRFdens_compute = True
        if depvars == 'WRFpos': WRFpos_compute = True
        if depvars == 'WRFtime': WRFtime_compute = True
        if depvars == 'WRFz': WRFz_compute = True
    else:
        depvars = diags[idiag].split('|')[1].split('@')
        if gen.searchInlist(depvars, 'WRFgeop'): WRFgeop_compute = True
        if gen.searchInlist(depvars, 'WRFp'): WRFp_compute = True
        if gen.searchInlist(depvars, 'WRFt'): WRFt_compute = True
        if gen.searchInlist(depvars, 'WRFrh'): WRFrh_compute = True
        if gen.searchInlist(depvars, 'WRFght'): WRFght_compute = True
        if gen.searchInlist(depvars, 'WRFdens'): WRFdens_compute = True
        if gen.searchInlist(depvars, 'WRFpos'): WRFpos_compute = True
        if gen.searchInlist(depvars, 'WRFtime'): WRFtime_compute = True
        if gen.searchInlist(depvars, 'WRFz'): WRFz_compute = True

# Dictionary with the new computed variables to be able to add them
dictcompvars = {}
if WRFgeop_compute:
    print '  ' + main + ': Retrieving geopotential value from WRF as PH + PHB'
    dimv = ncobj.variables['PH'].shape
    WRFgeop = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]

    # Attributes of the variable
    Vvals = gen.variables_values('WRFgeop')
    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}

if WRFp_compute:
    print '  ' + main + ': Retrieving pressure value from WRF as P + PB'
    dimv = ncobj.variables['P'].shape
    WRFp = ncobj.variables['P'][:] + ncobj.variables['PB'][:]

    # Attributes of the variable
    Vvals = gen.variables_values('WRFp')
    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}

if WRFght_compute:
    print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...' 
    WRFght = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]

    # Attributes of the variable
    Vvals = gen.variables_values('WRFght')
    dictcompvars['WRFgeop'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}

if WRFrh_compute:
    print '    ' + main + ": computing relative humidity from WRF as 'Tetens'" +     \
      ' equation (T,P) ...'
    p0=100000.
    p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
    tk = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
    qv = ncobj.variables['QVAPOR'][:]

    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)

    WRFrh = qv/data2

    # Attributes of the variable
    Vvals = gen.variables_values('WRFrh')
    dictcompvars['WRFrh'] = {'name': Vvals[0], 'standard_name': Vvals[1],            \
      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}

if WRFt_compute:
    print '    ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
    p0=100000.
    p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]

    WRFt = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)

    # Attributes of the variable
    Vvals = gen.variables_values('WRFt')
    dictcompvars['WRFt'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}

if WRFdens_compute:
    print '    ' + main + ': computing air density from WRF as ((MU + MUB) * ' +     \
      'DNW)/g ...'

# Just we need in in absolute values: Size of the central grid cell
##    dxval = ncobj.getncattr('DX')
##    dyval = ncobj.getncattr('DY')
##    mapfac = ncobj.variables['MAPFAC_M'][:]
##    area = dxval*dyval*mapfac

    mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
    dnw = ncobj.variables['DNW'][:]

    WRFdens = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]),        \
      dtype=np.float)
    levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)

    for it in range(mu.shape[0]):
        for iz in range(dnw.shape[1]):
            levval.fill(np.abs(dnw[it,iz]))
            WRFdens[it,iz,:,:] = levval
            WRFdens[it,iz,:,:] = mu[it,:,:]*WRFdens[it,iz,:,:]/grav

    # Attributes of the variable
    Vvals = gen.variables_values('WRFdens')
    dictcompvars['WRFdens'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}

if WRFpos_compute:
# WRF positions from the lowest-leftest corner of the matrix
    print '    ' + main + ': computing position from MAPFAC_M as sqrt(DY*j**2 + ' +  \
      'DX*x**2)*MAPFAC_M ...'

    mapfac = ncobj.variables['MAPFAC_M'][:]

    distx = np.float(ncobj.getncattr('DX'))
    disty = np.float(ncobj.getncattr('DY'))

    print 'distx:',distx,'disty:',disty

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

    WRFpos = np.zeros((dt, dy, dx), dtype=np.float)

    for i in range(1,dx):
        WRFpos[0,0,i] = distx*i/mapfac[0,0,i]
    for j in range(1,dy):
        i=0
        WRFpos[0,j,i] = WRFpos[0,j-1,i] + disty/mapfac[0,j,i]
        for i in range(1,dx):
#            WRFpos[0,j,i] = np.sqrt((disty*j)**2. + (distx*i)**2.)/mapfac[0,j,i]
#            WRFpos[0,j,i] = np.sqrt((disty*j)**2. + (distx*i)**2.)
             WRFpos[0,j,i] = WRFpos[0,j,i-1] + distx/mapfac[0,j,i]

    for it in range(1,dt):
        WRFpos[it,:,:] = WRFpos[0,:,:]

if WRFtime_compute:
    print '    ' + main + ': computing time from WRF as CFtime(Times) ...'

    refdate='19491201000000'
    tunitsval='minutes'

    timeobj = ncobj.variables['Times']
    timewrfv = timeobj[:]

    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 = timeobj.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

    # Attributes of the variable
    dictcompvars['WRFtime'] = {'name': 'time', 'standard_name': 'time',              \
      'long_name': 'time', 'units': tunits, 'calendar': 'gregorian'}

if WRFz_compute:
    print '  ' + main + ': Retrieving z: height above surface value from WRF as ' +  \
      'unstagger(PH + PHB)/9.8-hgt'
    dimv = ncobj.variables['PH'].shape
    WRFzg = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/9.8

    unzgd = (dimv[0], dimv[1]-1, dimv[2], dimv[3])
    unzg = np.zeros(unzgd, dtype=np.float)
    unzg = 0.5*(WRFzg[:,0:dimv[1]-1,:,:] + WRFzg[:,1:dimv[1],:,:])

    WRFz = np.zeros(unzgd, dtype=np.float)
    for iz in range(dimv[1]-1):
        WRFz[:,iz,:,:] = unzg[:,iz,:,:] - ncobj.variables['HGT'][:]

    # Attributes of the variable
    Vvals = gen.variables_values('WRFz')
    dictcompvars['WRFz'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}

### ## #
# Going for the diagnostics
### ## #
print '  ' + main + ' ...'
varsadd = []

for idiag in range(Ndiags):
    print '    diagnostic:',diags[idiag]
    diagn = diags[idiag].split('|')[0]
    depvars = diags[idiag].split('|')[1].split('@')
    if diags[idiag].split('|')[1].find('@') != -1:
        depvars = diags[idiag].split('|')[1].split('@')
        if depvars[0] == 'deaccum': diagn='deaccum'
        if depvars[0] == 'accum': diagn='accum'
        for depv in depvars:
            if not ncobj.variables.has_key(depv) and not                             \
              gen.searchInlist(NONcheckingvars, depv) and                            \
              not gen.searchInlist(methods, depv) and not depvars[0] == 'deaccum'    \
              and not depvars[0] == 'accum' and not depv[0:2] == 'z=':
                print errormsg
                print '  ' + main + ": file '" + opts.ncfile +                       \
                  "' does not have variable '" + depv + "' !!"
                quit(-1)
    else:
        depvars = diags[idiag].split('|')[1]
        if not ncobj.variables.has_key(depvars) and not                              \
          gen.searchInlist(NONcheckingvars, depvars) and                             \
          not gen.searchInlist(methods, depvars):
            print errormsg
            print '  ' + main + ": file '" + opts.ncfile +                           \
              "' does not have variable '" + depvars + "' !!"
            quit(-1)

    print "\n    Computing '" + diagn + "' from: ", depvars, '...'

# acraintot: accumulated total precipitation from WRF RAINC, RAINNC
    if diagn == 'ACRAINTOT':
            
        var0 = ncobj.variables[depvars[0]]
        var1 = ncobj.variables[depvars[1]]
        diagout = var0[:] + var1[:]

        dnamesvar = var0.dimensions
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnamesvar,nonvd): dvnamesvar.remove(nonvd)
            varsadd.append(nonvd)

        ncvar.insert_variable(ncobj, 'pracc', diagout, dnamesvar, dvnamesvar, newnc)

# accum: acumulation of any variable as (Variable, time [as [tunits] 
#   from/since ....], newvarname)
    elif diagn == 'accum':

        var0 = ncobj.variables[depvars[0]]
        var1 = ncobj.variables[depvars[1]]

        dnamesvar = var0.dimensions
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout, diagoutd, diagoutvd = diag.compute_accum(var0,dnamesvar,dvnamesvar)

        CFvarn = ncvar.variables_values(depvars[0])[0]

# Removing the flux
        if depvars[1] == 'XTIME':
            dtimeunits = var1.getncattr('description')
            tunits = dtimeunits.split(' ')[0]
        else:
            dtimeunits = var1.getncattr('units')
            tunits = dtimeunits.split(' ')[0]

        dtime = (var1[1] - var1[0])*timeunits_seconds(tunits)

        ncvar.insert_variable(ncobj, CFvarn + 'acc', diagout*dtime, diagoutd, diagoutvd, newnc)

# cllmh with cldfra, pres
    elif diagn == 'cllmh':
            
        var0 = ncobj.variables[depvars[0]]
        if depvars[1] == 'WRFp':
            var1 = WRFp
        else:
            var01 = ncobj.variables[depvars[1]]
            if len(size(var1.shape)) < len(size(var0.shape)):
                var1 = np.brodcast_arrays(var01,var0)[0]
            else:
                var1 = var01

        diagout, diagoutd, diagoutvd = diag.Forcompute_cllmh(var0,var1,dnames,dvnames)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        for nonvd in NONchkvardims:
            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)

        ncvar.insert_variable(ncobj, 'cll', diagout[0,:], diagoutd, diagoutvd, newnc)
        ncvar.insert_variable(ncobj, 'clm', diagout[1,:], diagoutd, diagoutvd, newnc)
        ncvar.insert_variable(ncobj, 'clh', diagout[2,:], diagoutd, diagoutvd, newnc)

# clt with cldfra
    elif diagn == 'clt':
            
        var0 = ncobj.variables[depvars]
        diagout, diagoutd, diagoutvd = diag.Forcompute_clt(var0,dnames,dvnames)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        for nonvd in NONchkvardims:
            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)
            
        ncvar.insert_variable(ncobj, 'clt', diagout, diagoutd, diagoutvd, newnc)

# deaccum: deacumulation of any variable as (Variable, time [as [tunits] 
#   from/since ....], newvarname)
    elif diagn == 'deaccum':

        var0 = ncobj.variables[depvars[1]]
        var1 = ncobj.variables[depvars[2]]

        dnamesvar = var0.dimensions
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout, diagoutd, diagoutvd = diag.compute_deaccum(var0,dnamesvar,dvnamesvar)

# Transforming to a flux
        if depvars[2] == 'XTIME':
            dtimeunits = var1.getncattr('description')
            tunits = dtimeunits.split(' ')[0]
        else:
            dtimeunits = var1.getncattr('units')
            tunits = dtimeunits.split(' ')[0]

        dtime = (var1[1] - var1[0])*timeunits_seconds(tunits)
        ncvar.insert_variable(ncobj, depvars[3], diagout/dtime, diagoutd, diagoutvd, newnc)

# LMDZrh (pres, t, r)
    elif diagn == 'LMDZrh':
            
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:]
        var2 = ncobj.variables[depvars[2]][:]

        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnames,dvnames)
        ncvar.insert_variable(ncobj, 'hur', diagout, diagoutd, diagoutvd, newnc)

# LMDZrhs (psol, t2m, q2m)
    elif diagn == 'LMDZrhs':
            
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:]
        var2 = ncobj.variables[depvars[2]][:]

        dnamesvar = ncobj.variables[depvars[0]].dimensions
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)

        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)

# mrso: total soil moisture SMOIS, DZS
    elif diagn == 'WRFmrso':
            
        var0 = ncobj.variables[depvars[0]][:]
        var10 = ncobj.variables[depvars[1]][:]
        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        var1 = var0.copy()*0.
        var2 = var0.copy()*0.+1.
        # Must be a better way....
        for j in range(var0.shape[2]):
          for i in range(var0.shape[3]):
              var1[:,:,j,i] = var10

        diagout, diagoutd, diagoutvd = diag.Forcompute_zint(var0, var1, var2,        \
          dnamesvar, dvnamesvar)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        diagoutvd = list(dvnames)
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)
        ncvar.insert_variable(ncobj, 'mrso', diagout, diagoutd, diagoutvd, newnc)

# mslp: mean sea level pressure (pres, psfc, terrain, temp, qv)
    elif diagn == 'mslp' or diagn == 'WRFmslp':
            
        var1 = ncobj.variables[depvars[1]][:]
        var2 = ncobj.variables[depvars[2]][:]
        var4 = ncobj.variables[depvars[4]][:]

        if diagn == 'WRFmslp':
            var0 = WRFp
            var3 = WRFt
            dnamesvar = ncobj.variables['P'].dimensions
        else:
            var0 = ncobj.variables[depvars[0]][:]
            var3 = ncobj.variables[depvars[3]][:]
            dnamesvar = ncobj.variables[depvars[0]].dimensions

        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout, diagoutd, diagoutvd = diag.compute_mslp(var0, var1, var2, var3, var4,    \
          dnamesvar, dvnamesvar)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        diagoutvd = list(dvnames)
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)
        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)

# OMEGAw (omega, p, t) from NCL formulation (https://www.ncl.ucar.edu/Document/Functions/Contributed/omega_to_w.shtml)
    elif diagn == 'OMEGAw':
            
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:]
        var2 = ncobj.variables[depvars[2]][:]

        dnamesvar = ncobj.variables[depvars[0]].dimensions
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout, diagoutd, diagoutvd = diag.compute_OMEGAw(var0,var1,var2,dnamesvar,dvnamesvar)

        ncvar.insert_variable(ncobj, 'wa', diagout, diagoutd, diagoutvd, newnc)

# raintot: instantaneous total precipitation from WRF as (RAINC + RAINC) / dTime
    elif diagn == 'RAINTOT':

        var0 = ncobj.variables[depvars[0]]
        var1 = ncobj.variables[depvars[1]]
        if depvars[2] != 'WRFtime':
            var2 = ncobj.variables[depvars[2]]
        else:
            var2 = np.arange(var0.shape[0], dtype=int)

        var = var0[:] + var1[:]

        dnamesvar = var0.dimensions
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout, diagoutd, diagoutvd = diag.compute_deaccum(var,dnamesvar,dvnamesvar)

# Transforming to a flux
        if var2.shape[0] > 1:
            if depvars[2] != 'WRFtime':
                dtimeunits = var2.getncattr('units')
                tunits = dtimeunits.split(' ')[0]
   
                dtime = (var2[1] - var2[0])*timeunits_seconds(tunits)
            else:
                var2 = ncobj.variables['Times']
                time1 = var2[0,:]
                time2 = var2[1,:]
                tmf1 = ''
                tmf2 = ''
                for ic in range(len(time1)):
                    tmf1 = tmf1 + time1[ic]
                    tmf2 = tmf2 + time2[ic]
                dtdate1 = dtime.datetime.strptime(tmf1,"%Y-%m-%d_%H:%M:%S")
                dtdate2 = dtime.datetime.strptime(tmf2,"%Y-%m-%d_%H:%M:%S")
                diffdate12 = dtdate2 - dtdate1
                dtime = diffdate12.total_seconds()
                print 'dtime:',dtime
        else:
            print warnmsg
            print '  ' + main + ": only 1 time-step for '" + diag + "' !!"
            print '    leaving a zero value!'
            diagout = var0[:]*0.
            dtime=1.

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        for nonvd in NONchkvardims:
            if gen.searchInlist(diagoutvd,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)
            
        ncvar.insert_variable(ncobj, 'pr', diagout/dtime, diagoutd, diagoutvd, newnc)

# rhs (psfc, t, q) from TimeSeries files
    elif diagn == 'TSrhs':
            
        p0=100000.
        var0 = ncobj.variables[depvars[0]][:]
        var1 = (ncobj.variables[depvars[1]][:])*(var0/p0)**(2./7.)
        var2 = ncobj.variables[depvars[2]][:]

        dnamesvar = ncobj.variables[depvars[0]].dimensions
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)

        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)

# slw: total soil liquid water SH2O, DZS
    elif diagn == 'WRFslw':
            
        var0 = ncobj.variables[depvars[0]][:]
        var10 = ncobj.variables[depvars[1]][:]
        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        var1 = var0.copy()*0.
        var2 = var0.copy()*0.+1.
        # Must be a better way....
        for j in range(var0.shape[2]):
          for i in range(var0.shape[3]):
              var1[:,:,j,i] = var10

        diagout, diagoutd, diagoutvd = diag.Forcompute_zint(var0, var1, var2,        \
          dnamesvar, dvnamesvar)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        diagoutvd = list(dvnames)
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)
        ncvar.insert_variable(ncobj, 'slw', diagout, diagoutd, diagoutvd, newnc)

# td (psfc, t, q) from TimeSeries files
    elif diagn == 'TStd' or diagn == 'td':
            
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:] - 273.15
        var2 = ncobj.variables[depvars[2]][:]

        dnamesvar = ncobj.variables[depvars[0]].dimensions
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)

        ncvar.insert_variable(ncobj, 'tds', diagout, diagoutd, diagoutvd, newnc)

# td (psfc, t, q) from TimeSeries files
    elif diagn == 'TStdC' or diagn == 'tdC':
            
        var0 = ncobj.variables[depvars[0]][:]
# Temperature is already in degrees Celsius
        var1 = ncobj.variables[depvars[1]][:]
        var2 = ncobj.variables[depvars[2]][:]

        dnamesvar = ncobj.variables[depvars[0]].dimensions
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout, diagoutd, diagoutvd = diag.compute_td(var0,var1,var2,dnamesvar,dvnamesvar)

        ncvar.insert_variable(ncobj, 'tds', diagout, diagoutd, diagoutvd, newnc)

# wds (u, v)
    elif diagn == 'TSwds' or diagn == 'wds' :
 
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:]

        dnamesvar = ncobj.variables[depvars[0]].dimensions
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout, diagoutd, diagoutvd = diag.compute_wds(var0,var1,dnamesvar,dvnamesvar)

        ncvar.insert_variable(ncobj, 'wds', diagout, diagoutd, diagoutvd, newnc)

# wss (u, v)
    elif diagn == 'TSwss' or diagn == 'wss':
            
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:]

        dnamesvar = ncobj.variables[depvars[0]].dimensions
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout, diagoutd, diagoutvd = diag.compute_wss(var0,var1,dnamesvar,dvnamesvar)

        ncvar.insert_variable(ncobj, 'wss', diagout, diagoutd, diagoutvd, newnc)

# turbulence (var)
    elif diagn == 'turbulence':

        var0 = ncobj.variables[depvars][:]

        dnamesvar = list(ncobj.variables[depvars].dimensions)
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout, diagoutd, diagoutvd = diag.compute_turbulence(var0,dnamesvar,dvnamesvar)
        valsvar = gen.variables_values(depvars)

        newvarn = depvars + 'turb'
        print main + '; Lluis newvarn:', newvarn
        ncvar.insert_variable(ncobj, newvarn, diagout, diagoutd, 
          diagoutvd, newnc)
        print main + '; Lluis variables:', newnc.variables.keys()
        varobj = newnc.variables[newvarn]
        attrv = varobj.long_name
        attr = varobj.delncattr('long_name')
        newattr = ncvar.set_attribute(varobj, 'long_name', attrv +                   \
          " Taylor decomposition turbulence term")

# WRFbils fom WRF as HFX + LH
    elif diagn == 'WRFbils':
            
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:]

        diagout = var0 + var1
        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        ncvar.insert_variable(ncobj, 'bils', diagout, dnamesvar, dvnamesvar, newnc)

# WRFcape_afwa CAPE, CIN, ZLFC, PLFC, LI following WRF 'phys/module_diaf_afwa.F' 
#   methodology as WRFt, WRFrh, WRFp, WRFgeop, HGT
    elif diagn == 'WRFcape_afwa':
        var0 = WRFt
        var1 = WRFrh
        var2 = WRFp
        dz = WRFgeop.shape[1]
        # de-staggering
        var3 = 0.5*(WRFgeop[:,0:dz-1,:,:]+WRFgeop[:,1:dz,:,:])
        var4 = ncobj.variables[depvars[4]][0,:,:]

        dnamesvar = list(ncobj.variables['T'].dimensions)
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout = np.zeros(var0.shape, dtype=np.float)
        diagout1, diagout2, diagout3, diagout4, diagout5, diagoutd, diagoutvd =      \
          diag.Forcompute_cape_afwa(var0, var1, var2, var3, var4, 3, dnamesvar,      \
          dvnamesvar)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)

        ncvar.insert_variable(ncobj, 'cape', diagout1, diagoutd, diagoutvd, newnc)
        ncvar.insert_variable(ncobj, 'cin', diagout2, diagoutd, diagoutvd, newnc)
        ncvar.insert_variable(ncobj, 'zlfc', diagout3, diagoutd, diagoutvd, newnc)
        ncvar.insert_variable(ncobj, 'plfc', diagout4, diagoutd, diagoutvd, newnc)
        ncvar.insert_variable(ncobj, 'li', diagout5, diagoutd, diagoutvd, newnc)

# WRFclivi WRF water vapour path WRFdens, QICE, QGRAUPEL, QHAIL
    elif diagn == 'WRFclivi':
            
        var0 = WRFdens
        qtot = ncobj.variables[depvars[1]]
        qtotv = qtot[:]
        Nspecies = len(depvars) - 2
        for iv in range(Nspecies):
            if ncobj.variables.has_key(depvars[iv+2]):
                var1 = ncobj.variables[depvars[iv+2]][:]
                qtotv = qtotv + var1

        dnamesvar = list(qtot.dimensions)
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout, diagoutd, diagoutvd = diag.compute_clivi(var0, qtotv, dnamesvar,dvnamesvar)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        diagoutvd = list(dvnames)
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)
        ncvar.insert_variable(ncobj, 'clivi', diagout, diagoutd, diagoutvd, newnc)

# WRFclwvl WRF water cloud-condensed path WRFdens, QCLOUD, QICE, QGRAUPEL, QHAIL
    elif diagn == 'WRFclwvl':
            
        var0 = WRFdens
        qtot = ncobj.variables[depvars[1]]
        qtotv = ncobj.variables[depvars[1]]
        Nspecies = len(depvars) - 2
        for iv in range(Nspecies):
            if ncobj.variables.has_key(depvars[iv+2]):
                var1 = ncobj.variables[depvars[iv+2]]
                qtotv = qtotv + var1[:]

        dnamesvar = list(qtot.dimensions)
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout, diagoutd, diagoutvd = diag.compute_clwvl(var0, qtotv, dnamesvar,dvnamesvar)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        diagoutvd = list(dvnames)
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)
        ncvar.insert_variable(ncobj, 'clwvl', diagout, diagoutd, diagoutvd, newnc)

# WRFgeop geopotential from WRF as PH + PHB
    elif diagn == 'WRFgeop':
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:]

        # de-staggering geopotential
        diagout0 = var0 + var1
        dt = diagout0.shape[0]
        dz = diagout0.shape[1]
        dy = diagout0.shape[2]
        dx = diagout0.shape[3]

        diagout = np.zeros((dt,dz-1,dy,dx), dtype=np.float)
        diagout = 0.5*(diagout0[:,1:dz,:,:]+diagout0[:,0:dz-1,:,:])

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        diagoutvd = list(dvnames)
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)

        ncvar.insert_variable(ncobj, 'zg', diagout, dnames, diagoutvd, newnc)

# WRFmslp_ptarget sea-level pressure following ECMWF method as PSFC, HGT, WRFt, WRFp, ZNU, ZNW
    elif diagn == 'WRFpsl_ecmwf':
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][0,:,:]
        var2 = WRFt[:,0,:,:]
        var4 = WRFp[:,0,:,:]
        var5 = ncobj.variables[depvars[4]][0,:]
        var6 = ncobj.variables[depvars[5]][0,:]

        # This is quite too appriximate!! passing pressure at half-levels to 2nd full 
        #   level, using eta values at full (ZNW) and half (ZNU) mass levels
        var3 = WRFp[:,0,:,:] + (var6[1] - var5[0])*(WRFp[:,1,:,:] - WRFp[:,0,:,:])/  \
          (var5[1]-var5[0])

        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout = np.zeros(var0.shape, dtype=np.float)
        diagout, diagoutd, diagoutvd = diag.Forcompute_psl_ecmwf(var0, var1, var2,   \
          var3, var4, dnamesvar, dvnamesvar)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)

        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)

# WRFmslp_ptarget sea-level pressure following ptarget method as WRFp, PSFC, WRFt, HGT, QVAPOR
    elif diagn == 'WRFpsl_ptarget':
        var0 = WRFp
        var1 = ncobj.variables[depvars[1]][:]
        var2 = WRFt
        var3 = ncobj.variables[depvars[3]][0,:,:]
        var4 = ncobj.variables[depvars[4]][:]

        dnamesvar = list(ncobj.variables[depvars[4]].dimensions)
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout = np.zeros(var0.shape, dtype=np.float)
        diagout, diagoutd, diagoutvd = diag.Forcompute_psl_ptarget(var0, var1, var2, \
          var3, var4, 700000., dnamesvar, dvnamesvar)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)

        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)

# WRFp pressure from WRF as P + PB
    elif diagn == 'WRFp':
            
        diagout = WRFp

        ncvar.insert_variable(ncobj, 'pres', diagout, dnames, dvnames, newnc)

# WRFpos 
    elif diagn == 'WRFpos':
            
        dnamesvar = ncobj.variables['MAPFAC_M'].dimensions
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        ncvar.insert_variable(ncobj, 'WRFpos', WRFpos, dnamesvar, dvnamesvar, newnc)

# WRFprw WRF water vapour path WRFdens, QVAPOR
    elif diagn == 'WRFprw':
            
        var0 = WRFdens
        var1 = ncobj.variables[depvars[1]]

        dnamesvar = list(var1.dimensions)
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout, diagoutd, diagoutvd = diag.compute_prw(var0, var1, dnamesvar,dvnamesvar)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        diagoutvd = list(dvnames)
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)
        ncvar.insert_variable(ncobj, 'prw', diagout, diagoutd, diagoutvd, newnc)

# WRFrh (P, T, QVAPOR)
    elif diagn == 'WRFrh':
            
        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        ncvar.insert_variable(ncobj, 'hur', WRFrh, dnames, dvnames, newnc)

# WRFrhs (PSFC, T2, Q2)
    elif diagn == 'WRFrhs':
            
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:]
        var2 = ncobj.variables[depvars[2]][:]

        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        diagout, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
        ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)

# rvors (u10, v10, WRFpos)
    elif diagn == 'WRFrvors':
            
        var0 = ncobj.variables[depvars[0]]
        var1 = ncobj.variables[depvars[1]]

        diagout = rotational_z(var0, var1, distx)

        dnamesvar = ncobj.variables[depvars[0]].dimensions
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        ncvar.insert_variable(ncobj, 'rvors', diagout, dnamesvar, dvnamesvar, newnc)

# WRFt (T, P, PB)
    elif diagn == 'WRFt':
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:]
        var2 = ncobj.variables[depvars[2]][:]

        p0=100000.
        p=var1 + var2

        WRFt = (var0 + 300.)*(p/p0)**(2./7.)

        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        diagoutvd = list(dvnames)
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)

        ncvar.insert_variable(ncobj, 'ta', WRFt, dnames, diagoutvd, newnc)

# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
    elif diagn == 'WRFua':
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:]
        var2 = ncobj.variables[depvars[2]][:]
        var3 = ncobj.variables[depvars[3]][:]

        # 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,dnames,dvnames)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        diagoutvd = list(dvnames)
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)

        ncvar.insert_variable(ncobj, 'ua', ua, dnames, diagoutvd, newnc)

# WRFua (U, V, SINALPHA, COSALPHA) to be rotated !!
    elif diagn == 'WRFva':
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:]
        var2 = ncobj.variables[depvars[2]][:]
        var3 = ncobj.variables[depvars[3]][:]

        # 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,dnames,dvnames)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        diagoutvd = list(dvnames)
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)
        ncvar.insert_variable(ncobj, 'va', va, dnames, diagoutvd, newnc)

# WRFtime
    elif diagn == 'WRFtime':
            
        diagout = WRFtime

        dnamesvar = ['Time']
        dvnamesvar = ['Times']

        ncvar.insert_variable(ncobj, 'time', diagout, dnamesvar, dvnamesvar, newnc)

# ws (U, V)
    elif diagn == 'ws':
            
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:]
        # 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],:])

        dnamesvar = ['Time','bottom_top','south_north','west_east']
        diagout = np.sqrt(unstgvar0*unstgvar0 + unstgvar1*unstgvar1)

#        dnamesvar = ncobj.variables[depvars[0]].dimensions
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        diagoutvd = list(dvnamesvar)
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnamesvar,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)
        ncvar.insert_variable(ncobj, 'ws', diagout, dnamesvar, diagoutvd, newnc)

# wss (u10, v10)
    elif diagn == 'wss':
            
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:]

        diagout = np.sqrt(var0*var0 + var1*var1)

        dnamesvar = ncobj.variables[depvars[0]].dimensions
        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)

        ncvar.insert_variable(ncobj, 'wss', diagout, dnamesvar, dvnamesvar, newnc)

# WRFheight height from WRF geopotential as WRFGeop/g
    elif diagn == 'WRFheight':
            
        diagout = WRFgeop/grav

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        diagoutvd = list(dvnames)
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)

        ncvar.insert_variable(ncobj, 'zhgt', diagout, dnames, diagoutvd, newnc)

# WRFheightrel relative-height from WRF geopotential as WRFgeop(PH + PHB)/g-HGT 'WRFheightrel|PH@PHB@HGT
    elif diagn == 'WRFheightrel':
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:]
        var2 = ncobj.variables[depvars[2]][:]

        dimz = var0.shape[1]
        diagout = np.zeros(tuple(var0.shape), dtype=np.float)
        for iz in range(dimz):
            diagout[:,iz,:,:] = (var0[:,iz,:,:]+ var1[:,iz,:,:])/grav - var2

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        diagoutvd = list(dvnames)
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)

        ncvar.insert_variable(ncobj, 'zhgtrel', diagout, dnames, diagoutvd, newnc)

# WRFzmla_gen generic boundary layer hieght computation from WRF theta, QVAPOR, WRFgeop, HGT, 
    elif diagn == 'WRFzmlagen':
        var0 = ncobj.variables[depvars[0]][:]+300.
        var1 = ncobj.variables[depvars[1]][:]
        dimz = var0.shape[1]
        var2 = WRFgeop[:,1:dimz+1,:,:]/9.8
        var3 = ncobj.variables[depvars[3]][0,:,:]

        diagout, diagoutd, diagoutvd = diag.Forcompute_zmla_gen(var0,var1,var2,var3, \
          dnames,dvnames)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)

        ncvar.insert_variable(ncobj, 'zmla', diagout, diagoutd, diagoutvd, newnc)

# WRFzwind wind extrapolation at a given height using power law computation from WRF 
#   U, V, WRFz, U10, V10, SINALPHA, COSALPHA, z=[zval]
    elif diagn == 'WRFzwind':
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:]
        var2 = WRFz
        var3 = ncobj.variables[depvars[3]][:]
        var4 = ncobj.variables[depvars[4]][:]
        var5 = ncobj.variables[depvars[5]][0,:,:]
        var6 = ncobj.variables[depvars[6]][0,:,:]
        var7 = np.float(depvars[7].split('=')[1])

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

        diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwind(unvar0,      \
          unvar1, var2, var3, var4, var5, var6, var7, dnames, dvnames)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)

        ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
        ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)

# WRFzwind wind extrapolation at a given hieght using logarithmic law computation 
#   from WRF U, V, WRFz, U10, V10, SINALPHA, COSALPHA, z=[zval]
    elif diagn == 'WRFzwind_log':
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:]
        var2 = WRFz
        var3 = ncobj.variables[depvars[3]][:]
        var4 = ncobj.variables[depvars[4]][:]
        var5 = ncobj.variables[depvars[5]][0,:,:]
        var6 = ncobj.variables[depvars[6]][0,:,:]
        var7 = np.float(depvars[7].split('=')[1])

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

        diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwind_log(unvar0,  \
          unvar1, var2, var3, var4, var5, var6, var7, dnames, dvnames)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)

        ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
        ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)

# WRFzwindMO wind extrapolation at a given height computation using Monin-Obukhow 
#   theory from WRF UST, ZNT, RMOL, U10, V10, SINALPHA, COSALPHA, z=[zval]
#   NOTE: only useful for [zval] < 80. m
    elif diagn == 'WRFzwindMO':
        var0 = ncobj.variables[depvars[0]][:]
        var1 = ncobj.variables[depvars[1]][:]
        var2 = ncobj.variables[depvars[2]][:]
        var3 = ncobj.variables[depvars[3]][:]
        var4 = ncobj.variables[depvars[4]][:]
        var5 = ncobj.variables[depvars[5]][0,:,:]
        var6 = ncobj.variables[depvars[6]][0,:,:]
        var7 = np.float(depvars[7].split('=')[1])

        diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwindMO(var0, var1,\
          var2, var3, var4, var5, var6, var7, dnames, dvnames)

        # Removing the nonChecking variable-dimensions from the initial list
        varsadd = []
        for nonvd in NONchkvardims:
            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
            varsadd.append(nonvd)

        ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
        ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)

    else:
        print errormsg
        print '  ' + main + ": diagnostic '" + diagn + "' not ready!!!"
        print '    available diagnostics: ', availdiags
        quit(-1)

    newnc.sync()
    # Adding that additional variables required to compute some diagnostics which
    #   where not in the original file
    for vadd in varsadd:
        if not gen.searchInlist(newnc.variables.keys(),vadd):
            attrs = dictcompvars[vadd]
            vvn = attrs['name']
            if not gen.searchInlist(newnc.variables.keys(), vvn):
                iidvn = dvnames.index(vadd)
                dnn = dnames[iidvn]
                if vadd == 'WRFtime':
                    dvarvals = WRFtime[:]
                newvar = newnc.createVariable(vvn, 'f8', (dnn))
                newvar[:] = dvarvals
                for attn in attrs.keys():
                    if attn != 'name':
                        attv = attrs[attn]
                        ncvar.set_attribute(newvar, attn, attv)

#   end of diagnostics

# Global attributes
##
ncvar.add_global_PyNCplot(newnc, main, None, '2.0')

gorigattrs = ncobj.ncattrs()
for attr in gorigattrs:
    attrv = ncobj.getncattr(attr)
    atvar = ncvar.set_attribute(newnc, attr, attrv)

ncobj.close()
newnc.close()

print '\n' + main + ': successfull writting of diagnostics file "' + ofile + '" !!!'