Changeset 1687 in lmdz_wrf

Timestamp:

Dec 5, 2017, 1:11:15 PM (7 years ago)

Author:

lfita

Message:

Adding:

• WRF computinf of Individual wind components
• `C_diagnostic': Class to compute generic variables
• `W_diagnostic': Class to compute WRF diagnostics variables

File:

• trunk/tools/diag_tools.py (modified) (7 diffs)

trunk/tools/diag_tools.py

@@ -16 +16 @@
 # 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_wds: Function to compute the wind direction
 # compute_wss: Function to compute the wind speed
 # 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
 # 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
+# W_diagnostic: Class to compute WRF diagnostics variables
 
 # Others just providing variable values
@@ -33 +39 @@
 # 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
@@ -52 +60 @@
 # Constants
 grav = fdef.module_definitions.grav
+
+# Available WRFiag
+Wavailablediags = ['p', 'ta', 'td', 'ua', 'va', 'uas', 'vas', 'wd', 'ws', 'zg']
+
+# Available General diagnostics
+Cavailablediags = ['td', 'wd', 'ws']
 
 # Gneral information
@@ -968 +982 @@
     return times
 
+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
@@ -1003 +1083 @@
     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
@@ -1072 +1206 @@
     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
+
 def var_td(t, p, qv):
     """ Function to compute dew-point air temperature from temperature and pressure values
@@ -1092 +1262 @@
 
     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 == 'td':
+            """ Computing dew-point temperature 
+            """
+            vn = 'td'
+            CF3Dvars = ['ta', 'plev', 'hur']
+            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']][:]
+            hur = ncobj.variables[vars3D['hur']][:]
+
+            self.values = var_td(ta, p, hur)
+            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 == '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'][:]
+            hur = ncobj.variables['QVAPOR'][:]
+    
+            vals, dims, vdims = compute_WRFtd(ta, p, hur, 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)