Index: /palm/trunk/UTIL/WRF_interface/README.txt
===================================================================
--- /palm/trunk/UTIL/WRF_interface/README.txt	(revision 4766)
+++ /palm/trunk/UTIL/WRF_interface/README.txt	(revision 4766)
@@ -0,0 +1,116 @@
+Scripts for processing of WRF and CAMx files to PALM dynamic driver.
+Version: v.3.1, 20201102
+
+Usage: palm_dynamic -c <config_name> [-w]
+
+The script requires name of the case configuration on the command line.
+The corresponding configuration files are placed in subdirectory
+"configuration" and they are named <config_name>.conf. Examples of
+the config files are supplied with the scripts, all configuration
+options with are in the file palm_dynamic_defaults.py as well as
+at the end of this readme file. The values which agree with defaults
+need not be present in the user config. The file palm_dynamic_init.py
+contains setting and calculation of standard initialization values
+for particular system and can be adjusted. The optional parameter -w
+allows to skip horizontal and vertical interpolation in case it is
+already done.
+
+The scripts were implemented with regards to good portability and they
+only depend on standard and well-known Python libraries.
+Needed modules are:
+
+- numpy   (https://pypi.org/project/numpy)
+- scipy   (https://pypi.org/project/scipy)
+- pyproj  (https://pypi.org/project/pyproj)
+- netCDF4 (https://pypi.org/project/netCDF4)
+- metpy   (https://unidata.github.io/MetPy)
+
+First four modules are included in all major linux distributions. In other cases,
+they can be installed via pip (e.g. "pip3 install netCDF4").
+
+The MetPy module is a well-known and widely used project for meteorological
+calculations in Python maintained by Unidata. It can also be installed using
+pip ("pip3 install MetPy").
+
+In the current version, the only supported projection in WRF is Lambert
+conformal conic, which is WRF default and recommended projection for
+mid-latitudes. The CAMx output files, which are used optionally for the
+chemical initial and boundary conditions, are expected to be the result of
+a coupled CAMx simulation using the same horizontal grid as WRF, although this
+is not strictly required.
+
+The scripts support both variants of WRF vertical levels - the sigma levels
+(default until WRF version 3.*) and the hybrid levels (default since WRF 4.*).
+However, it is necessary to correctly configure this option via the setting
+"wrf_hybrid_levs = True/False".
+
+CONFIGURATION
+
+Configuration files have syntax:
+<parameter1> = <value1>
+<parameter2> = <value2>
+...
+Parameters can be string, integer, float, logical or list of these values.
+The <value> terms have Python syntax, the config file needs to contain only
+parameters which differ from default.
+
+Description of the particular configuration options are (defaults are in parenthesis):
+
+# 1. Domain and case related config
+domain              name of the simulation case ("")
+resolution          name of the particular domain resolution scenario ("")
+scenario            name of the individual scenario in the case ("")
+nested_domain       False indicates parent and True nested domain. (False)
+dynamic_driver_file file name of output dynamic driver (""). The default value
+                    "" means the standard name assigned in the palm_dynamic_init.
+grid_from_static    value True - the grid parameters are imported from the static
+                    driver, False - they are prescribed in the config (True)
+static_driver_file  file name of the static driver in case of grid_from_static ("").
+                    The default value "" means the standard name assigned in init.
+proj_palm           reference coordinate system of PALM simulation ("EPSG:32633")
+proj_wgs84          reference coordinate system of lon-lat projection ("EPSG:4326")
+dz                  height of the PALM vertical grid layer (0.0). The default
+                    value dz = 0.0 means dz is assigned from dx.
+nz                  number of vertical layers of PALM domain (200)
+dz_stretch_level    height in meters from which stretching of vertical levels
+                    starts in PALM (5000.0)
+dz_stretch_factor   coefficient of the stretching of the vertical layers in PALM (1.0)
+dz_max              max height of the stretched vertical layers (100.0)
+origin_time         origin time of the PALM simulation in the format
+                    YYYY-MM-DD hh:mm:ss (""). The default value "" means that
+                    the value is read from the global attribute of the static driver.
+simulation_hours    extent of the simulation in hours (24)
+
+# 2. WRF and CAMx related configurations
+wrf_hybrid_levs     True means hybrid levels in WRF files, False means sigma levels (True).
+vinterp_terrain_smoothing
+                    the standard deviation for Gaussian kernel of smoothing method
+                    of the PALM terrain for WRF vertical interpolation to avoid sharp
+                    horizontal gradients. Value None disables the smoothing. (None)
+wrf_dir_name        files path of the wrf and camx input files (""). The default value ""
+                    means that the standard path will be calculated in the init.
+wrf_file_mask       file mask of the wrf input files  ("wrfout_*.e000")
+radiation_from_wrf  enable or disable processing of radiation from WRF files (True).
+wrf_rad_file_mask   file mask of the wrf radiation input files ("auxhist6_*").
+                    The default setting reads radiation from WRF auxiliary
+                    history files. This setting allows to use finer time step for WRF
+                    radiation outputs than for other values.
+radiation_smoothing_distance
+                    smoothing distance for radiation values in m (10000.0).
+camx_file_mask      file mask of the CAMx input files ("CAMx.*.nc").
+species_names       PALM names of the chemical species (['NO','NO2','O3','PM10','PM25']).
+                    The default value can be used for phstatp and phstatp2 chemical
+                    mechanisms in PALM. The mapping and recalculation from the CAMx
+                    species to PALM species is defined in the variable camx_conversions
+                    (see palm_dynamic_init.py). The current mapping is adjusted for
+                    CAMx v6.50 with CB05 + PM (CF,SOAP2.1,ISORROPIA) chemistry mechanisms.
+                    In the case other combinations, it may be necessary to adjust
+                    this mapping.
+
+# 3. horizontal parameters of the PALM domain which have to be set in case
+#    of grid_from_static = False
+nx, ny              number of horizontal grids of the domain in x and y directions
+dx, dy              grid cell size of the domain in x and y directions
+origin_x, origin_y  origin x and y of the domain
+origin_z            origin of the domain in the vertical direction
+
Index: /palm/trunk/UTIL/WRF_interface/configurations/dejvice_validation_summer_02.conf
===================================================================
--- /palm/trunk/UTIL/WRF_interface/configurations/dejvice_validation_summer_02.conf	(revision 4766)
+++ /palm/trunk/UTIL/WRF_interface/configurations/dejvice_validation_summer_02.conf	(revision 4766)
@@ -0,0 +1,17 @@
+# Simulation for summer validation in Prague-Dejvice
+# see GMD: https://doi.org/10.5194/gmd-2020-175
+# nested domain with 2 m resolution
+
+domain = 'dejvice'
+resolution = '02'
+scenario = 'validation_summer'
+
+# vertical grid
+nz = 120  # z in grids
+nested_domain = True
+
+origin_time = "2018-07-19 00:00:00"
+simulation_hours = 72
+
+# wrf and camx files
+wrf_dir_name = '/storage/data/palm/dejvice/validation/wrf/summer'
Index: /palm/trunk/UTIL/WRF_interface/configurations/dejvice_validation_summer_10.conf
===================================================================
--- /palm/trunk/UTIL/WRF_interface/configurations/dejvice_validation_summer_10.conf	(revision 4766)
+++ /palm/trunk/UTIL/WRF_interface/configurations/dejvice_validation_summer_10.conf	(revision 4766)
@@ -0,0 +1,19 @@
+# Simulation for summer validation in Prague-Dejvice
+# see GMD: https://doi.org/10.5194/gmd-2020-175
+# parent domain with 10 m resolution
+
+domain = 'dejvice'
+resolution = '10'
+scenario = 'validation_summer'
+
+# vertical grid
+nz = 160  # z in grids
+dz_stretch_level = 250.0 # in meters
+dz_stretch_factor = 1.08
+dz_max = 20.0
+
+origin_time = "2018-07-19 00:00:00"
+simulation_hours = 72
+
+# wrf and camx files
+wrf_dir_name = '/storage/data/palm/dejvice/validation/wrf/summer'
Index: /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic.py
===================================================================
--- /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic.py	(revision 4766)
+++ /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic.py	(revision 4766)
@@ -0,0 +1,528 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+#------------------------------------------------------------------------------#
+#
+# Scripts for processing of WRF and CAMx files to PALM dynamic driver
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+#
+# Copyright 2018-2020 Institute of Computer Science
+#                     of the Czech Academy of Sciences, Prague
+# Authors: Krystof Eben, Jaroslav Resler, Pavel Krc
+#
+#------------------------------------------------------------------------------#
+'''Scripts for processing of WRF and CAMx files to PALM dynamic driver.
+
+Usage: palm_dynamic -c <config_name> [-w]
+
+Script requires name of the configuration on command line.
+The corresponding case configuration files are placed in subdirectory
+"configuration" and the are named <config_name>.conf. The template
+of the config file is in the file palm_dynamic_defaults.py, the values
+which agree with defaults need not be present in the user config.
+The file palm_dynamic_init.py contains setting and calculation of
+standard initialization values for particular system and can be adjusted.
+The optional parameter -w allows to skip horizontal and vertical
+interpolation in case it is already done.
+'''
+__version__ = '3.1'
+
+import sys
+import getopt
+import os.path
+from datetime import datetime, timedelta
+import numpy as np
+from pyproj import Proj, transform
+import glob
+
+import netCDF4
+
+import palm_dynamic_config
+
+##################################
+# read configuration from command
+# line and parse it
+def print_help():
+    print()
+    print(__doc__)
+    print('Version: '+__version__)
+
+# set commandline variables
+configname = ''
+wrf_interpolation = True
+if sys.argv[0].endswith('pydevconsole.py'):
+    # we are in Pycharm developlent console - set testing data
+    configname = 'evropska_12s_basecase'
+else:
+    try:
+        opts, args = getopt.getopt(sys.argv[1:],"hc:w)",["help=","config=","write="])
+        #print(opts)
+    except getopt.GetoptError as err:
+        print("Error:", err)
+        print_help()
+        sys.exit(2)
+    # parse options
+    for opt, arg in opts:
+        if opt in ("-h", "--help"):
+            print_help()
+            sys.exit(0)
+        elif opt in ("-c", "--config"):
+            configname = arg
+        elif opt in ("-w", "--write"):
+            wrf_interpolation = False
+
+if configname == '':
+    # script requires config name
+    print("This script requires input parameter -c <config_name>")
+    print_help()
+    exit(2)
+
+# Load config defaults, config file and standard init into the config module
+palm_dynamic_config.configure(configname)
+# Import values (including loaded) from config module into globals
+from palm_dynamic_config import *
+
+# Other modules are imported *AFTER* the config module has been filled with
+# values (so that they can correctly import * from palm_dynamic_config).
+import palm_wrf_utils
+from palm_dynamic_output import palm_dynamic_output
+import palm_dynamic_camx
+
+
+print('Domain name: ', domain)
+print('Resolution name: ', resolution)
+print('Scenario name: ', scenario)
+print('Read domain from static:', grid_from_static)
+if grid_from_static:
+    print('STATIC driver input file: ', static_driver_file)
+print('WRF and CAMx file path:', wrf_dir_name)
+print('WRF file mask:', wrf_file_mask)
+print('CAMx file mask:', camx_file_mask)
+print('Simulation start time:', origin_time)
+print('Simulation hours:', simulation_hours)
+
+
+###########################
+# domain parameters
+
+if grid_from_static:
+    # get parameters of the horizontal domain from PALM STATIC driver
+    try:
+        ncs = netCDF4.Dataset(static_driver_file, "r", format="NETCDF4")
+    except Exception as err:
+        print("Error opening static driver file:")
+        print(static_driver_file)
+        print(err)
+        sys.exit(1)
+    # get horizontal structure of the domain
+    nx = ncs.dimensions['x'].size
+    ny = ncs.dimensions['y'].size
+    dx = ncs.variables['x'][:][1] - ncs.variables['x'][:][0]
+    dy = ncs.variables['y'][:][1] - ncs.variables['y'][:][0]
+    origin_x = ncs.getncattr('origin_x')
+    origin_y = ncs.getncattr('origin_y')
+    origin_z = ncs.getncattr('origin_z')
+    if origin_time == '':
+        # origin_time is not prvided in configuration - read it from static driver
+        origin_time = ncs.getncattr('origin_time')
+
+    # create vertical structure of the domain
+    # calculate dz in case dz is not supplied (dz<=0)
+    if dz <= 0.0:
+        print("dz = 0.0: set dz = dx")
+        dz = dx
+    if nz <= 0:
+        print("nz > 0 needs to be set in config")
+        sys.exit(1)
+
+    # read terrain height (relative to origin_z) and
+    # calculate and check the height of the surface canopy layer
+    if 'zt' in ncs.variables.keys():
+        terrain_rel = ncs.variables['zt'][:]
+    else:
+        terrain_rel = np.zeros([ny,nx])
+    th = np.ceil(terrain_rel / dz)
+    # building height
+    if 'buildings_3d' in ncs.variables.keys():
+        print(np.argmax(a != 0, axis=0))
+        bh3 = ncs.variables['buildings_3d'][:]
+        # minimum index of nonzeo value along inverted z
+        bh = np.argmax(bh3[::-1], axis=0)
+        # inversion back and masking grids with no buildings
+        bh = bh3.shape[0] - bh
+        bh[np.max(bh3, axis=0) == 0] = 0
+    elif 'buildings_2d' in ncs.variables.keys():
+        bh = ncs.variables['buildings_2d'][:]
+        bh[bh.mask] = 0
+        bh = np.ceil(bh / dz)
+    else:
+        bh = np.zeros([ny,nx])
+    # plant canopy height
+    if 'lad' in ncs.variables.keys():
+        lad3 = ncs.variables['lad'][:]
+        # replace non-zero values with 1
+        lad3[lad3 != 0] = 1
+        # minimum index of nonzeo value along inverted z
+        lad = np.argmax(lad3[::-1], axis=0)
+        # inversion back and masking grids with no buildings
+        lad = lad3.shape[0] - lad
+        lad[np.max(lad3, axis=0) == 0] = 0
+    else:
+        lad = np.zeros([ny,nx])
+    # calculate maximum of surface canopy layer
+    nscl = max(np.amax(th+bh),np.amax(th+lad))
+    # check nz with ncl
+    if nz < nscl + nscl_free:
+        print('nz has to be higher than ', nscl + nscl_free)
+        print('nz=', nz, ', dz=', dz, ', number of scl=', nscl, ',nscl_free=', nscl_free)
+    if dz_stretch_level <= (nscl + nscl_free) * dz:
+        print('stretching has to start in higher level than ', (nscl + nscl_free) * dz)
+        print('dz_stretch_level=', dz_stretch_level, ', nscl=', nscl, ', nscl_free=', nscl_free, ', dz=', dz)
+    if 'soil_moisture_adjust' in ncs.variables.keys():
+        soil_moisture_adjust = ncs.variables['soil_moisture_adjust'][:]
+    else:
+        soil_moisture_adjust = np.ones(shape=(ny, nx), dtype=float)
+    # close static driver nc file
+    ncs.close()
+else:
+    #TODO check all necessary parameters are set and are correct
+    #
+    # initialize necessary arrays
+    try:
+        terrain = np.zeros(shape=(ny, nx), dtype=float)
+        soil_moisture_adjust = np.ones(shape=(ny, nx), dtype=float)
+    except:
+        print('domain parameters nx, ny have to be set in configuration')
+        sys.exit(1)
+
+# absolute terrain needed for vertical interpolation of wrf data
+terrain = terrain_rel + origin_z
+
+# print domain parameters and check ist existence in caso of setup from config
+try:
+    print('Domain parameters:')
+    print('nx, ny, nz:', nx, ny, nz)
+    print('dx, dy, dz:', dx, dy, dz)
+    print('origin_x, origin_y:', origin_x, origin_y)
+    print('Base of domain is in level origin_z:', origin_z)
+except:
+    print('domain parameters have to be read from static driver or set in configuration')
+    sys.exit(1)
+
+# centre of the domain (needed for ug,vg calculation)
+xcent = origin_x + nx * dx / 2.0
+ycent = origin_y + ny * dy / 2.0
+# WGS84 projection for transformation to lat-lon
+inproj = Proj('+init='+proj_palm)
+print('inproj', inproj)
+lonlatproj = Proj('+init='+proj_wgs84)
+print('lonlatproj', lonlatproj)
+cent_lon, cent_lat = transform(inproj, lonlatproj, xcent, ycent)
+print('xcent, ycent:',xcent, ycent)
+print('cent_lon, cent_lat:', cent_lon, cent_lat)
+# prepare target grid
+irange = origin_x + dx * (np.arange(nx, dtype='f8') + .5)
+jrange = origin_y + dy * (np.arange(ny, dtype='f8') + .5)
+palm_grid_y, palm_grid_x = np.meshgrid(jrange, irange, indexing='ij')
+palm_grid_lon, palm_grid_lat = transform(inproj, lonlatproj, palm_grid_x, palm_grid_y)
+
+######################################
+# build structure of vertical layers
+# remark:
+# PALM input requires nz=ztop in PALM
+# but the output file in PALM has max z higher than z in PARIN.
+# The highest levels in PALM are wrongly initialized !!!
+#####################################
+# fill out z_levels
+z_levels = np.zeros(nz,dtype=float)
+z_levels_stag = np.zeros(nz-1,dtype=float)
+dzs = dz
+z_levels[0] = dzs/2.0
+for i in range(nz-1):
+    z_levels[i+1] = z_levels[i] + dzs
+    z_levels_stag[i] = (z_levels[i+1]+z_levels[i])/2.0
+    dzso = dzs
+    if z_levels[i+1] + dzs >= dz_stretch_level:
+        dzs = min(dzs * dz_stretch_factor, dz_max)
+ztop = z_levels[-1] + dzs / 2.
+print('z:',z_levels)
+print('zw:',z_levels_stag)
+
+######################################
+# get time extent of the PALM simulation
+#####################################
+# get complete list of wrf files
+wrf_file_list = glob.glob(os.path.join(wrf_dir_name, wrf_file_mask))
+# get simulation origin and final time as datetime
+start_time = datetime.strptime(origin_time, '%Y-%m-%d %H:%M:%S')
+end_time = start_time + timedelta(hours=simulation_hours)
+end_time_rad = end_time
+print('PALM simulation extent', start_time, end_time, simulation_hours)
+if nested_domain:
+    print('Nested domain - process only initialization.')
+    print('Set end_time = start_time')
+    end_time = start_time
+
+# get wrf times and sort wrf files by time
+print('Analyse WRF files dates:')
+file_times = []
+for wrf_file in wrf_file_list:
+    # get real time from wrf file
+    nc_wrf = netCDF4.Dataset(wrf_file, "r", format="NETCDF4")
+    ta = nc_wrf.variables['Times'][:]
+    t = ta.tobytes().decode("utf-8")
+    td = datetime.strptime(t, '%Y-%m-%d_%H:%M:%S')
+    print(os.path.basename(wrf_file), ': ', td)
+    file_times.append((td,wrf_file))
+file_times.sort()
+times = []
+wrf_files = []
+for tf in file_times:
+    if end_time is None or tf[0] <= end_time:
+        times.append(tf[0])
+        wrf_files.append(tf[1])
+#print('PALM output times:')
+#print('\n'.join('{}'.format(t) for t in times))
+print('PALM output times:', ', '.join('{}'.format(t) for t in times))
+
+if not times.__contains__(start_time):
+    print('WRF files does not contain PALM origin_time timestep - cannot process!')
+    exit(1)
+
+if not times.__contains__(end_time):
+    print('WRF files does not contain PALM end_time timestep - cannot process!')
+    exit(1)
+
+start_index = times.index(start_time)
+end_index = times.index(end_time)
+
+if not nested_domain and end_index-start_index < simulation_hours:
+    print('Number of WRF files does not aggre with number of simulation hours')
+    exit(1)
+
+print('PALM simulation timestep number are from ', start_index, ' to ', end_index, ' WRF file.')
+
+# create list of processed wrf files
+wrf_files_proc = []
+for wf in wrf_files[start_index:end_index+1]:
+    wrf_files_proc.append(wf)
+# id of process files
+simul_id = domain + "_d" + resolution
+
+######################################
+#VERTICAL AND HORIZONTAL INTERPOLATION
+######################################
+# get hydro and soil variables contained in wrf files
+shvars = []
+z_soil_levels = []
+nc_wrf = netCDF4.Dataset(wrf_files_proc[0], "r", format="NETCDF4")
+try:
+    # hydro variables
+    shvars = sorted(set(['QVAPOR', 'QCLOUD', 'QRAIN', 'QICE', 'QSNOW', 'QGRAUP']).intersection(nc_wrf.variables.keys()))
+    # soil layers
+    if 'ZS' in nc_wrf.variables.keys():
+        z_soil_levels = nc_wrf.variables['ZS'][0,:].data.tolist()
+finally:
+    nc_wrf.close()
+print('Hydro variables in wrf files:', '+'.join(shvars))
+
+print('Start of vertical and horizontal interpolation of inputs to the PALM domain.')
+interp_files = []
+regridder = None
+for wrf_file in wrf_files_proc:
+    print ("Input wrf file: ", wrf_file)
+    #
+    hinterp_file = wrf_file+"_"+simul_id+'.hinterp'
+    hinterp_log = wrf_file + "_" + simul_id + '.hinterp.log'
+    vinterp_file = wrf_file+"_"+simul_id+'.interp'
+    interp_files.append(vinterp_file)
+    if wrf_interpolation:
+        try:
+             os.remove(vinterp_file)
+             os.remove(hinterp_file)
+             os.remove(hinterp_log)
+        except:
+              pass
+
+        print('Horizontal interpolation...')
+        f_wrf = netCDF4.Dataset(wrf_file, 'r')
+        try:
+            if not regridder:
+                trans_wrf = palm_wrf_utils.WRFCoordTransform(f_wrf)
+                palm_in_wrf_y, palm_in_wrf_x = trans_wrf.latlon_to_ji(
+                                                palm_grid_lat, palm_grid_lon)
+                regridder = palm_wrf_utils.BilinearRegridder(
+                        palm_in_wrf_x, palm_in_wrf_y, preloaded=True)
+                regridder_u = palm_wrf_utils.BilinearRegridder(
+                        palm_in_wrf_x+.5, palm_in_wrf_y, preloaded=True)
+                regridder_v = palm_wrf_utils.BilinearRegridder(
+                        palm_in_wrf_x, palm_in_wrf_y+.5, preloaded=True)
+
+            f_out = netCDF4.Dataset(hinterp_file, 'w', format='NETCDF4')
+            try:
+                # dimensions
+                f_out.createDimension('Time', None)
+                for d in 'bottom_top bottom_top_stag soil_layers_stag'.split():
+                    f_out.createDimension(d, len(f_wrf.dimensions[d]))
+                f_out.createDimension('west_east', len(irange))
+                f_out.createDimension('south_north', len(jrange))
+
+                # copied vars
+                for varname in 'PH PHB HGT T W TSLB SMOIS MU MUB P PB PSFC'.split():
+                    v_wrf = f_wrf.variables[varname]
+                    v_out = f_out.createVariable(varname, 'f4', v_wrf.dimensions)
+                    v_out[:] = regridder.regrid(v_wrf[...,regridder.ys,regridder.xs])
+
+                # U and V have special treatment (unstaggering)
+                v_out = f_out.createVariable('U', 'f4', ('Time', 'bottom_top', 'south_north', 'west_east'))
+                v_out[:] = regridder_u.regrid(f_wrf.variables['U'][...,regridder_u.ys,regridder_u.xs])
+                v_out = f_out.createVariable('V', 'f4', ('Time', 'bottom_top', 'south_north', 'west_east'))
+                v_out[:] = regridder_v.regrid(f_wrf.variables['V'][...,regridder_v.ys,regridder_v.xs])
+
+                # calculated SPECHUM
+                v_out = f_out.createVariable('SPECHUM', 'f4', ('Time', 'bottom_top', 'south_north', 'west_east'))
+                vdata = regridder.regrid(f_wrf.variables[shvars[0]][...,regridder.ys,regridder.xs])
+                for vname in shvars[1:]:
+                    vdata += regridder.regrid(f_wrf.variables[vname][...,regridder.ys,regridder.xs])
+                v_out[:] = vdata
+            finally:
+                f_out.close()
+        finally:
+            f_wrf.close()
+
+        print('Vertical interpolation...')
+        palm_wrf_utils.palm_wrf_vertical_interp(hinterp_file, vinterp_file, wrf_file, z_levels,
+                                                z_levels_stag, z_soil_levels, origin_z, terrain,
+                                                wrf_hybrid_levs, vinterp_terrain_smoothing)
+
+if radiation_from_wrf:
+    print('Start processing of radiation inputs from the WRF radiation files.')
+    # get available times from wrf radiation output files and sort files by time
+    # get complete list of wrf radiation files
+    rad_file_list = glob.glob(os.path.join(wrf_dir_name, wrf_rad_file_mask))
+    rad_file_times = []
+    print('Analyse WRF radiation files dates:')
+    for rad_file in rad_file_list:
+        # get real time from wrf radiation file
+        nc_rad = netCDF4.Dataset(rad_file, "r", format="NETCDF4")
+        ta = nc_rad.variables['Times'][:]
+        t = ta.tobytes().decode("utf-8")
+        td = datetime.strptime(t, '%Y-%m-%d_%H:%M:%S')
+        print(os.path.basename(rad_file), ': ', td)
+        rad_file_times.append((td, rad_file))
+    rad_file_times.sort()
+    rad_times = []
+    rad_files = []
+    for tf in rad_file_times:
+        if end_time_rad is None or tf[0] <= end_time_rad:
+            rad_times.append(tf[0])
+            rad_files.append(tf[1])
+    # check list of available times
+    if not rad_times.__contains__(start_time):
+        print('WRF radiation files does not contain PALM origin_time timestep - cannot process!')
+        exit(1)
+    if not rad_times.__contains__(end_time_rad):
+        print('WRF radiation files does not contain PALM end_time_rad timestep - cannot process!')
+        exit(1)
+    rad_start_index = rad_times.index(start_time)
+    rad_end_index = rad_times.index(end_time_rad)
+    print('PALM radiation timestep numbers are from ', rad_start_index, ' to ', rad_end_index, ' WRF radiation file.')
+    # get radiation timestep
+    rad_timestep = (rad_times[rad_start_index+1] - rad_times[rad_start_index]).total_seconds()
+    print('Timestep of the wrf radiation data is ', rad_timestep,' seconds.')
+    # check consistency of the series
+    for i in range(rad_start_index, rad_end_index):
+        if (rad_times[i+1] - rad_times[i]).total_seconds() != rad_timestep:
+            print('Inconsistent timeline of radiation inputs! Cannot process.')
+            print(rad_times[rad_start_index:rad_end_index])
+            exit(1)
+    # create list of processed wrf radiation files and times
+    rad_files_proc = []
+    rad_times_proc = []
+    for i in range(rad_start_index, rad_end_index + 1):
+        rad_files_proc.append(rad_files[i])
+        rad_times_proc.append((rad_times[i]-rad_times[rad_start_index]).total_seconds())
+    print("Radiation times are ", rad_times_proc)
+    # process radiation inputs
+    rad_swdown = []
+    rad_lwdown = []
+    rad_swdiff = []
+    nfiles = 0
+    for rad_file in rad_files_proc:
+        print ("Input wrf radiation file: ", rad_file)
+        ncf = netCDF4.Dataset(rad_file, "r", format="NETCDF4")
+        try:
+            if nfiles == 0:
+                # create list of (i,j) indices used for smoothing of the radiation
+                print('Build list of indices for radiation smoothig.')
+                palmproj = Proj(init=proj_palm)
+                lonlatproj = Proj(init=proj_wgs84)
+                rad_ind = []
+                for i in range(ncf.dimensions['west_east'].size):
+                    for j in range(ncf.dimensions['south_north'].size):
+                        # transfrom lat-lon to PALM domain in meters
+                        lonij = ncf.variables['XLONG'][0,j,i]
+                        latij = ncf.variables['XLAT'][0,j,i]
+                        xij, yij = transform(lonlatproj, palmproj, lonij, latij)
+                        if abs(xij-xcent)<=radiation_smoothing_distance and \
+                            abs(yij-ycent)<=radiation_smoothing_distance:
+                            rad_ind.append([i,j])
+                ngrids = len(rad_ind)
+            # process wrf radiation file
+            nfiles += 1
+            swdown = 0.0
+            lwd = 0.0
+            swddif = 0.0
+            for ij in rad_ind:
+                swdown = swdown + ncf.variables['SWDOWN'][0,ij[1],ij[0]]
+                lwd    = lwd    + ncf.variables['GLW'][0,ij[1],ij[0]]
+                swddif = swddif + ncf.variables['SWDDIF'][0,ij[1],ij[0]]
+            rad_swdown.append(swdown / ngrids)
+            rad_lwdown.append(lwd / ngrids)
+            rad_swdiff.append(swddif / ngrids)
+        finally:
+            ncf.close()
+    # create list of all radiation values
+    rad_values_proc = [rad_swdown, rad_lwdown, rad_swdiff]
+else:
+    rad_times_proc = []
+    rad_values_proc = []
+
+camx_interp_fname = None
+if not nested_domain:
+    # process camx files
+    camx_file_list = glob.glob(os.path.join(wrf_dir_name, camx_file_mask))
+    if camx_file_list:
+        print('Processing CAMx input files: {0}'.format(', '.join(camx_file_list)))
+        camx_interp_fname = os.path.join(wrf_dir_name, 'CAMx_{0}_interp.nc'.format(simul_id))
+
+        palm_dynamic_camx.process_files(camx_file_list, camx_interp_fname,
+                palm_grid_lat, palm_grid_lon, terrain_rel, z_levels,
+                times[start_index:end_index+1], species_names,
+                camx_conversions, camx_helpers)
+
+# ===================== CREATE NETCDF DRIVER ==============================
+# calculate relative times from simulation start
+times_sec = []
+for t in times[start_index:end_index+1]:
+    times_sec.append((t-start_time).total_seconds())
+# collect dimension sizes
+dimensions = {'zdim': nz, 'zwdim': nz-1, 'zsoildim': len(z_soil_levels), 'xdim': nx, 'xudim': nx-1, 'ydim': ny, 'yvdim': ny-1}
+# process interpolated files to dynamic driver
+palm_dynamic_output(wrf_files, interp_files, camx_interp_fname, dynamic_driver_file, times_sec, dimensions,
+                    z_levels, z_levels_stag, ztop, z_soil_levels, dx, dy, cent_lon, cent_lat,
+                    rad_times_proc, rad_values_proc, soil_moisture_adjust, nested_domain)
+
+print('Creation of dynamic driver finished.')
Index: /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic_camx.py
===================================================================
--- /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic_camx.py	(revision 4766)
+++ /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic_camx.py	(revision 4766)
@@ -0,0 +1,195 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+#------------------------------------------------------------------------------#
+#
+# Scripts for processing of WRF and CAMx files to PALM dynamic driver
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+#
+# Copyright 2018-2020 Institute of Computer Science
+#                     of the Czech Academy of Sciences, Prague
+# Authors: Krystof Eben, Jaroslav Resler, Pavel Krc
+#
+#------------------------------------------------------------------------------#
+'''
+Load and process initial and boundary conditions from CAMx NetCDF files
+'''
+
+import sys
+import datetime
+import re
+import numpy as np
+import netCDF4
+from metpy.interpolate import interpolate_1d
+
+import palm_wrf_utils
+
+_na = np.newaxis
+re_num = re.compile(r'[0-9\.]+')
+
+class HelperNotFound(Exception):
+    def __init__(self, name):
+        self.name = name
+
+class Helpers(object):
+    '''Pre-loaded helper variables from CAMx'''
+    def __getattr__(self, name):
+        try:
+            return self.__dict__[name]
+        except KeyError:
+            raise HelperNotFound(name)
+
+def tflag(data, req_dts):
+    assert len(data.shape) == 3 and data.shape[2] == 2
+    xdate = data[:,0,0]
+    xtime = data[:,0,1]
+
+    # Verify that dates are equal for each variable
+    assert (data[:,:,0] == xdate[:,_na]).all()
+    assert (data[:,:,1] == xtime[:,_na]).all()
+
+    dts = []
+    for i in range(len(xdate)):
+        dt = datetime.datetime.strptime(
+            '{0:07d} {1:06d}'.format(xdate[i], xtime[i]), '%Y%j %H%M%S')
+        try:
+            ireq = req_dts[dt]
+        except KeyError:
+            continue
+        dts.append((ireq, i))
+    return dts
+
+def load_conversion(spc, f, sl, hlp, camx_vars, camx_units, formula, **kwargs):
+    loaded_vars = []
+    for varname, unit in zip(camx_vars, camx_units):
+        try:
+            v = f.variables[varname]
+        except KeyError:
+            print('Skipping {0}, because input variable {1}({2}) is missing.'.format(
+                spc, varname, unit))
+            return None
+        if getattr(v, 'units', None) != unit:
+            print('Skipping {0}, because input variable {1} has wrong unit ({2} <> {3}).'.format(
+                spc, varname, getattr(v, 'units', None), unit))
+            return None
+
+        print('Loading variable {0}.'.format(varname))
+        loaded_vars.append(v[sl])
+    try:
+        val = formula(*(loaded_vars + [hlp]))
+    except HelperNotFound as e:
+        print('Skipping {0} - missing helper variable {1}.'.format(spc, e.name))
+        return None
+
+    return val
+
+def process_tstep(f, itf, regridder, lay_height, fout, itout, z_levels,
+        vars_remaining, filled, conversions, helpers):
+
+    # Load helper vars for this timestep
+    hlp = Helpers()
+    for helper_name, helper in helpers:
+        data = load_conversion(helper_name, f,
+                (itf,slice(None),regridder.ys,regridder.xs), hlp, **helper)
+        if data is not None:
+            setattr(hlp, helper_name, data)
+
+    # Load all usable vars for this timestep, regrid horizontally
+    varmeta = []
+    vardata = []
+    for spc in list(vars_remaining):
+        conv = conversions[spc]
+        data = load_conversion(spc, f, (itf,slice(None),regridder.ys,regridder.xs),
+                hlp, **conv)
+        if data is None:
+            continue
+
+        data = regridder.regrid(data)
+        vardata.append(np.r_[data[0:1], data]) #add peg below
+        varmeta.append((spc, conv['output_unit']))
+
+    # Perform vertical interpolation on all currently loaded vars at once
+    print('Interpolating vertically...')
+    vinterp = interpolate_1d(z_levels, lay_height, *vardata)
+    if len(vardata) == 1:
+        # return_list_always=True argument is only in later versions of MetPy
+        vinterp = [vinterp]
+    del vardata
+    for (vn, vu), vd in zip(varmeta, vinterp):
+        v = fout.variables[vn]
+        v[itout] = vd
+        v.units = vu
+        filled[vn][itout] = True
+
+def process_files(camx_file_list, camx_interp_fname, palm_grid_lat,
+        palm_grid_lon, terrain_rel, z_levels, times, species_names,
+        conversions, helpers):
+
+    terrain_shift = terrain_rel[_na,:,:]
+    lowest_layer = np.zeros(((1,) + palm_grid_lat.shape), dtype='f4')
+    lowest_layer[:] = -999.
+    tindex = dict((dt, i) for i, dt in enumerate(times))
+    filled = {}
+
+    with netCDF4.Dataset(camx_interp_fname, 'w', format='NETCDF4') as fout:
+        fout.createDimension('time', len(times))
+        fout.createDimension('z', len(z_levels))
+        fout.createDimension('y', palm_grid_lat.shape[0])
+        fout.createDimension('x', palm_grid_lat.shape[1])
+        for vn in species_names:
+            fout.createVariable(vn, 'f4', ('time', 'z', 'y', 'x'))
+            filled[vn] = [False] * len(times)
+
+        for fname in sorted(camx_file_list):
+            with netCDF4.Dataset(fname) as f:
+                dts = tflag(f.variables['TFLAG'][:], tindex)
+                if dts:
+                    print('Processing CAMx file {0}.'.format(fname))
+
+                    # preprare projection
+                    trans = palm_wrf_utils.CAMxCoordTransform(f)
+                    palm_in_camx_y, palm_in_camx_x = trans.latlon_to_ji(
+                                                    palm_grid_lat, palm_grid_lon)
+                    regridder = palm_wrf_utils.BilinearRegridder(
+                            palm_in_camx_x, palm_in_camx_y, preloaded=True)
+
+                    # locate layer heights
+                    try:
+                        vz = f.variables['z']
+                    except KeyError:
+                        print('Loading heights from separate file')
+                        with open(fname+'.heights') as fh:
+                            fix_hgt = np.array(list(map(float, re_num.findall(fh.read())))) * 1000. #orig in km
+                            fix_hgt = fix_hgt[:,_na,_na]
+                    else:
+                        print('Loading heights from variable z')
+                        fix_hgt = None
+
+                    for itout, itf in dts:
+                        print('Timestep {0}'.format(itout))
+                        vars_remaining = [vn for vn, vf in filled.items() if not vf[itout]]
+
+                        lay_height = fix_hgt if fix_hgt is not None else (
+                                regridder.regrid(vz[itf,:,regridder.ys,regridder.xs]))
+                        lay_height = np.r_[lowest_layer, lay_height + terrain_shift]
+                                                #add 1 pegging layer always below
+                        process_tstep(f, itf, regridder, lay_height, fout, itout, z_levels,
+                                vars_remaining, filled, conversions, helpers)
+                else:
+                    print('Skipping CAMx file {0} - no required times.'.format(fname))
+
+    if not all(all(vf) for vf in filled.values()):
+        sys.exit('CAMx data not complete - missing some variables/timesteps: {0}'
+                .format(filled))
Index: /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic_config.py
===================================================================
--- /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic_config.py	(revision 4766)
+++ /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic_config.py	(revision 4766)
@@ -0,0 +1,65 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+#------------------------------------------------------------------------------#
+#
+# Scripts for processing of WRF and CAMx files to PALM dynamic driver
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+#
+# Copyright 2018-2020 Institute of Computer Science
+#                     of the Czech Academy of Sciences, Prague
+# Authors: Krystof Eben, Jaroslav Resler, Pavel Krc
+#
+#------------------------------------------------------------------------------#
+'''Configuration module.
+
+Configuration options are sourced into this module's globals.
+'''
+
+import os.path
+from pathlib import Path
+import inspect
+
+
+# Just for PyCharm and similar IDEs to allow autocompletion from config values
+if False:
+    from palm_dynamic_defaults import *
+    from palm_dynamic_init import *
+
+def configure(configname):
+    global dir_scripts
+    # get path of the palm_dynamic script to source default config and init
+    dir_scripts = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
+    print('Running palm_dynamic from:', dir_scripts)
+    # use config defaults
+    configdefaultsfile = os.path.join(dir_scripts, "palm_dynamic_defaults.py")
+    print('Default case config: ', configdefaultsfile)
+    # user config file is located in configurations directory
+    configfile = os.path.join(os.path.abspath(Path(dir_scripts).parent), "configurations", configname + '.conf')
+    print('User case config:', configfile)
+    # initialization of standard parameters done in script
+    standardinitfile = os.path.join(dir_scripts, "palm_dynamic_init.py")
+    print('Standard initialization: ', standardinitfile)
+    # check existence of the supplied config file
+    if not os.path.isfile(configfile):
+        print("Config file " + configfile + " does not exists!")
+        print_help()
+        exit(2)
+    # read default config values
+    exec(open(configdefaultsfile).read(), globals())
+    # read user configuration
+    exec(open(configfile).read(), globals())
+    # perform the standard initialization
+    exec(open(standardinitfile).read(), globals())
Index: /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic_defaults.py
===================================================================
--- /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic_defaults.py	(revision 4766)
+++ /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic_defaults.py	(revision 4766)
@@ -0,0 +1,55 @@
+
+# config template with config defaults
+# not a separate module, only sourced from palm_dynamic_config
+
+# PALM case, domain, and configuration parameters
+# (only an example, needs to be rewriten in user config)
+domain = ''
+resolution = ''
+scenario = ''
+nested_domain = False
+
+# file name of output dynamic driver ("" means the standard name)
+dynamic_driver_file = ""
+
+# import grid parameters for dynamic driver from static driver
+grid_from_static = True
+# file name of static driver ("" means the standard name)
+static_driver_file = ""
+# reference coordinate system of PALM simulation
+proj_palm = "EPSG:32633"
+# projection lon-lat
+proj_wgs84 = 'EPSG:4326'
+
+# vertical grid
+# layer height (dz = 0.0 means dz is assigned from dx)
+dz = 0.0
+# we need description of the PALM vertical structure for dynamic driver
+nz = 200  # z in grids
+dz_stretch_level = 5000.0 # in meters
+dz_stretch_factor = 1.0
+dz_max = 100.0
+
+# time origin and extent of the simulation (format YYYY-MM-DD hh:mm:ss)
+origin_time = ""
+simulation_hours = 24
+
+# WRF related configurations
+wrf_hybrid_levs = True
+# Smoothing of PALM terrain for WRF vertical interpolation to avoid sharp
+# horizontal gradients. None = off
+vinterp_terrain_smoothing = None
+
+# wrf and camx input files path and default file mask
+wrf_dir_name = ""  # "" means that standard path will be calculated in the init
+wrf_file_mask = "wrfout_*.e000"
+# process radiation from wrf
+radiation_from_wrf = True
+wrf_rad_file_mask = "auxhist6_*"
+# smoothing distance for radiation
+radiation_smoothing_distance = 10000.0
+
+# chemical initial and boundary conditions - CAMx nesting
+camx_file_mask = 'CAMx.*.nc'
+# PALM species names - version for phstatp2 mechanism
+species_names = ['NO', 'NO2', 'O3', 'PM10', 'PM25']
Index: /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic_init.py
===================================================================
--- /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic_init.py	(revision 4766)
+++ /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic_init.py	(revision 4766)
@@ -0,0 +1,122 @@
+#!/usr/bin/python3
+# -*- coding: utf-8 -*-
+
+#------------------------------------------------------------------------------#
+#
+# Scripts for processing of WRF and CAMx files to PALM dynamic driver
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+#
+# Copyright 2018-2020 Institute of Computer Science
+#                     of the Czech Academy of Sciences, Prague
+# Authors: Krystof Eben, Jaroslav Resler, Pavel Krc
+#
+#------------------------------------------------------------------------------#
+'''
+This file provides initialization of the basic variables and structures base on
+the default and user config values. It can be adjusted to the needs of particular
+user system or structure of the data storage.
+'''
+
+import os
+from pathlib import Path
+import numpy as np
+
+# paths of directories
+dir_base = os.path.abspath(Path(dir_scripts).parent.parent)
+print('Base dir:', dir_base)
+dir_in = os.path.join(dir_base, domain)
+print('Input dir:', dir_in)
+if scenario == '':
+    dir_in_scen = dir_in
+else:
+    dir_in_scen = os.path.join(dir_in,scenario)
+    print('Scenario input dir:', dir_in_scen)
+dir_out = os.path.join(dir_base, domain)
+print('Output dir:', dir_out)
+
+# file names of PALM PIDS drivers
+# extension of file name for scenario
+scenario_ext = ("_"+scenario if scenario != "" else "")
+# static driver netcdf file name
+if static_driver_file == "":
+    static_driver_file = os.path.join(dir_out, domain+"_static_driver"+"_d"+resolution+scenario_ext+".nc")
+# dynamic driver netcdf file name
+if dynamic_driver_file == "":
+    dynamic_driver_file = os.path.join(dir_out, domain+"_dynamic_driver"+"_d"+resolution+scenario_ext+".nc")
+
+# parameters of dynamic driver
+# minimal number of free surface canopy layers above top of terrain with building and plant canopy
+nscl_free = 3
+# default path of wrf files in case it is not set in user config
+if wrf_dir_name == '':
+    wrf_dir_name = os.path.join(dir_in, 'wrf')
+
+# Settings for geostrophic wind
+gw_gfs_margin_deg = 5. #smoothing area in degrees lat/lon
+gw_wrf_margin_km = 10. #smoothing area in km
+#gw_alpha = .143 #GW vertical interpolation by power law
+gw_alpha = 1. #ignore wind power law, interpolate linearly
+
+
+# chemical initial and boundary coords - CAMx nesting
+# CAMx variable conversions
+camx_conversions = {
+        'NO':  dict(
+            camx_vars = ['NO'],
+            camx_units = ['ppmv'],
+            output_unit = 'ppm',
+            formula = lambda no, hlp: no,
+            ),
+        'NO2':  dict(
+            camx_vars = ['NO2'],
+            camx_units = ['ppmv'],
+            output_unit = 'ppm',
+            formula = lambda no2, hlp: no2,
+            ),
+        'NOX': dict(
+            camx_vars=['NO', 'NO2'],
+            camx_units=['ppmv', 'ppmv'],
+            output_unit='ppm',
+            formula=lambda no, no2, hlp: no + no2,
+        ),
+        'O3':  dict(
+            camx_vars = ['O3'],
+            camx_units = ['ppmv'],
+            output_unit = 'ppm',
+            formula = lambda o3, hlp: o3,
+            ),
+        'PM10':  dict(
+            camx_vars = ['CPRM'],
+            camx_units = ['micrograms m-3'],
+            output_unit = 'kg/m3',
+            formula = lambda cprm, hlp: (cprm+hlp.pm25) * 1e-9,
+            ),
+        'PM25':  dict(
+            camx_vars = [],
+            camx_units = [],
+            output_unit = 'kg/m3',
+            formula = lambda hlp: hlp.pm25 * 1e-9,
+            ),
+        }
+
+camx_helpers = [
+        ('pm25',  dict(
+            camx_vars = 'PSO4 PNO3 PNH4 POA PEC FPRM SOA1 SOA2 SOA3 SOA4 SOPA SOPB'.split(),
+            camx_units = ['micrograms m-3']*12,
+            formula = lambda *args: np.sum(args[:-1], axis=0), #last arg is hlp
+            )),
+        ]
+
+
Index: /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic_output.py
===================================================================
--- /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic_output.py	(revision 4766)
+++ /palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic_output.py	(revision 4766)
@@ -0,0 +1,378 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+#------------------------------------------------------------------------------#
+#
+# Scripts for processing of WRF and CAMx files to PALM dynamic driver
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+#
+# Copyright 2018-2020 Institute of Computer Science
+#                     of the Czech Academy of Sciences, Prague
+# Authors: Krystof Eben, Jaroslav Resler, Pavel Krc
+#
+#------------------------------------------------------------------------------#
+'''
+This file creates and writes the dynamic driver netcdf file based on preprepared
+transformed and interpolated wrf and camx files.
+'''
+
+import os
+import time
+import numpy as np
+import netCDF4
+from palm_wrf_utils import palm_wrf_gw
+
+
+def palm_dynamic_output(wrf_files, interp_files, camx_interp_fname, dynamic_driver_file, times_sec,
+                        dimensions, z_levels, z_levels_stag, ztop,
+                        z_soil_levels, dx, dy, lon_center, lat_center,
+                        rad_times_proc, rad_values_proc, sma, nested_domain):
+
+    print('Processing interpolated files to dynamic driver')
+    # dimension of the time coordinate
+    dimtimes = len(times_sec)
+    # other coordinates
+    dimnames = ['z', 'zw', 'zsoil', 'x','xu', 'y', 'yv']                # z  height agl in m, zw staggered, zsoil 4 lev from wrf 
+    dimsize_names = ['zdim' , 'zwdim', 'zsoildim', 'xdim', 'xudim', 'ydim', 'yvdim']    # palm: zw = z - 1
+    x = np.arange(dx, dimensions['xdim']*dx+dx, dx)     # clean this
+    print('dimension of x:', len(x))
+    y = np.arange(dy, dimensions['ydim']*dy+dy, dy)     # clean this
+    print('dimension of y:', len(y))
+    # fill values
+    fillvalue_float = float(-9999.0)
+    # check out file and remove for sure
+    try:
+        os.remove(dynamic_driver_file)
+    except:
+        pass
+    # create netcdf out file
+    print('Driver file:', dynamic_driver_file)
+    outfile = netCDF4.Dataset(dynamic_driver_file, "w", format="NETCDF4" )
+    try:
+        # time dimension and variable
+        outfile.createDimension('time', dimtimes)
+        _val_times =  outfile.createVariable('time',"f4", ("time"))
+        # other dimensions and corresponding variables
+        for _dim in zip(dimnames, dimsize_names):  #range (len(dimnames)):
+            print(_dim[0],_dim[1], dimensions[_dim[1]])
+            outfile.createDimension(_dim[0], dimensions[_dim[1]])
+        _val_z_levels = outfile.createVariable('z',"f4", ("z"))
+        _val_z_levels_stag = outfile.createVariable('zw',"f4", ("zw"))
+        _val_z_soil_levels = outfile.createVariable('zsoil',"f4", ("zsoil"))
+        _val_y = outfile.createVariable('y',"f4", ("y"))
+        _val_x = outfile.createVariable('x',"f4", ("x"))
+
+        # prepare influx/outflux area sizes
+        zstag_all = np.r_[0., z_levels_stag, ztop]
+        zwidths = zstag_all[1:] - zstag_all[:-1]
+        print('zwidths', zwidths)
+        areas_xb = np.zeros((len(z_levels), 1))
+        areas_xb[:,0] = zwidths * dy
+        areas_yb = np.zeros((len(z_levels), 1))
+        areas_yb[:,0] = zwidths * dx
+        areas_zb = dx*dy
+        area_boundaries = (areas_xb.sum()*dimensions['ydim']*2
+                + areas_yb.sum()*dimensions['xdim']*2
+                + areas_zb*dimensions['xdim']*dimensions['ydim'])
+
+        # write values for coordinates
+        _val_times[:] = times_sec[:]
+        _val_z_levels[:] = z_levels[:]
+        _val_z_levels_stag[:] = z_levels_stag[:]
+        _val_z_soil_levels[:] = z_soil_levels[:]
+        _val_y[:] = y[:]
+        _val_x[:] = x[:]
+
+        # initialization of the variables and setting of init_* variables
+        print("Processing initialization from file", interp_files[0])
+        # open corresponding
+        infile = netCDF4.Dataset(interp_files[0], "r", format="NETCDF4")
+        try:
+            # open variables in the input file
+            init_atmosphere_pt = infile.variables['init_atmosphere_pt']
+            init_atmosphere_qv = infile.variables['init_atmosphere_qv']
+            init_atmosphere_u = infile.variables['init_atmosphere_u']
+            init_atmosphere_v = infile.variables['init_atmosphere_v']
+            init_atmosphere_w = infile.variables['init_atmosphere_w']
+            init_soil_m = infile.variables['init_soil_m']
+            init_soil_t = infile.variables['init_soil_t']
+
+            # create netcdf structure
+            _val_init_atmosphere_pt = outfile.createVariable('init_atmosphere_pt', "f4", ("z", "y", "x"),
+                                                             fill_value=fillvalue_float)
+            _val_init_atmosphere_pt.setncattr('lod', 2)
+            _val_init_atmosphere_qv = outfile.createVariable('init_atmosphere_qv', "f4", ("z", "y", "x"),
+                                                             fill_value=fillvalue_float)
+            _val_init_atmosphere_qv.setncattr('lod', 2)
+            _val_init_atmosphere_u = outfile.createVariable('init_atmosphere_u', "f4", ("z", "y", "xu"),
+                                                            fill_value=fillvalue_float)
+            _val_init_atmosphere_u.setncattr('lod', 2)
+            _val_init_atmosphere_v = outfile.createVariable('init_atmosphere_v', "f4", ("z", "yv", "x"),
+                                                            fill_value=fillvalue_float)
+            _val_init_atmosphere_v.setncattr('lod', 2)
+            _val_init_atmosphere_w = outfile.createVariable('init_atmosphere_w', "f4", ("zw", "y", "x"),
+                                                            fill_value=fillvalue_float)
+            _val_init_atmosphere_w.setncattr('lod', 2)
+            _val_init_soil_t = outfile.createVariable('init_soil_t', "f4", ("zsoil", "y", "x"), fill_value=fillvalue_float)
+            _val_init_soil_t.setncattr('lod', 2)
+            _val_init_soil_m = outfile.createVariable('init_soil_m', "f4", ("zsoil", "y", "x"), fill_value=fillvalue_float)
+            _val_init_soil_m.setncattr('lod', 2)
+            # time dependent variables
+            if not nested_domain:
+                # SURFACE PRESSURE
+                _val_surface_forcing_surface_pressure = outfile.createVariable('surface_forcing_surface_pressure', "f4",
+                                                                               ("time"))
+                # BOUNDARY - vertical slices from left, right, south, north, top
+                varname = 'pt'
+                _val_ls_forcing_pt_left = outfile.createVariable('ls_forcing_left_' + varname, "f4", ("time", "z", "y"),
+                                                                 fill_value=fillvalue_float)
+                _val_ls_forcing_pt_right = outfile.createVariable('ls_forcing_right_' + varname, "f4", ("time", "z", "y"),
+                                                                  fill_value=fillvalue_float)
+                _val_ls_forcing_pt_south = outfile.createVariable('ls_forcing_south_' + varname, "f4", ("time", "z", "x"),
+                                                                  fill_value=fillvalue_float)
+                _val_ls_forcing_pt_north = outfile.createVariable('ls_forcing_north_' + varname, "f4", ("time", "z", "x"),
+                                                                  fill_value=fillvalue_float)
+                _val_ls_forcing_pt_top = outfile.createVariable('ls_forcing_top_' + varname, "f4", ("time", "y", "x"),
+                                                                fill_value=fillvalue_float)
+
+                varname = 'qv'
+                _val_ls_forcing_qv_left = outfile.createVariable('ls_forcing_left_' + varname, "f4", ("time", "z", "y"),
+                                                                 fill_value=fillvalue_float)
+                _val_ls_forcing_qv_right = outfile.createVariable('ls_forcing_right_' + varname, "f4", ("time", "z", "y"),
+                                                                  fill_value=fillvalue_float)
+                _val_ls_forcing_qv_south = outfile.createVariable('ls_forcing_south_' + varname, "f4", ("time", "z", "x"),
+                                                                  fill_value=fillvalue_float)
+                _val_ls_forcing_qv_north = outfile.createVariable('ls_forcing_north_' + varname, "f4", ("time", "z", "x"),
+                                                                  fill_value=fillvalue_float)
+                _val_ls_forcing_qv_top = outfile.createVariable('ls_forcing_top_' + varname, "f4", ("time", "y", "x"),
+                                                                fill_value=fillvalue_float)
+
+                varname = 'u'
+                _val_ls_forcing_u_left = outfile.createVariable('ls_forcing_left_' + varname, "f4", ("time", "z", "y"),
+                                                                fill_value=fillvalue_float)
+                _val_ls_forcing_u_right = outfile.createVariable('ls_forcing_right_' + varname, "f4", ("time", "z", "y"),
+                                                                 fill_value=fillvalue_float)
+                _val_ls_forcing_u_south = outfile.createVariable('ls_forcing_south_' + varname, "f4", ("time", "z", "xu"),
+                                                                 fill_value=fillvalue_float)
+                _val_ls_forcing_u_north = outfile.createVariable('ls_forcing_north_' + varname, "f4", ("time", "z", "xu"),
+                                                                 fill_value=fillvalue_float)
+                _val_ls_forcing_u_top = outfile.createVariable('ls_forcing_top_' + varname, "f4", ("time", "y", "xu"),
+                                                               fill_value=fillvalue_float)
+
+                varname = 'v'
+                _val_ls_forcing_v_left = outfile.createVariable('ls_forcing_left_' + varname, "f4", ("time", "z", "yv"),
+                                                                fill_value=fillvalue_float)
+                _val_ls_forcing_v_right = outfile.createVariable('ls_forcing_right_' + varname, "f4", ("time", "z", "yv"),
+                                                                 fill_value=fillvalue_float)
+                _val_ls_forcing_v_south = outfile.createVariable('ls_forcing_south_' + varname, "f4", ("time", "z", "x"),
+                                                                 fill_value=fillvalue_float)
+                _val_ls_forcing_v_north = outfile.createVariable('ls_forcing_north_' + varname, "f4", ("time", "z", "x"),
+                                                                 fill_value=fillvalue_float)
+                _val_ls_forcing_v_top = outfile.createVariable('ls_forcing_top_' + varname, "f4", ("time", "yv", "x"),
+                                                               fill_value=fillvalue_float)
+
+                varname = 'w'
+                _val_ls_forcing_w_left = outfile.createVariable('ls_forcing_left_' + varname, "f4", ("time", "zw", "y"),
+                                                                fill_value=fillvalue_float)
+                _val_ls_forcing_w_right = outfile.createVariable('ls_forcing_right_' + varname, "f4", ("time", "zw", "y"),
+                                                                 fill_value=fillvalue_float)
+                _val_ls_forcing_w_south = outfile.createVariable('ls_forcing_south_' + varname, "f4", ("time", "zw", "x"),
+                                                                 fill_value=fillvalue_float)
+                _val_ls_forcing_w_north = outfile.createVariable('ls_forcing_north_' + varname, "f4", ("time", "zw", "x"),
+                                                                 fill_value=fillvalue_float)
+                _val_ls_forcing_w_top = outfile.createVariable('ls_forcing_top_' + varname, "f4", ("time", "y", "x"),
+                                                               fill_value=fillvalue_float)
+
+                # geostrophic wind
+                _val_ls_forcing_ug = outfile.createVariable('ls_forcing_ug', "f4", ("time", "z"),
+                                                            fill_value=fillvalue_float)
+                _val_ls_forcing_vg = outfile.createVariable('ls_forcing_vg', "f4", ("time", "z"),
+                                                            fill_value=fillvalue_float)
+                time.sleep(1)
+
+            # write values for initialization variables
+            _val_init_atmosphere_pt[:, :, :] = init_atmosphere_pt[0, :, :, :]
+            _val_init_atmosphere_qv[:, :, :] = init_atmosphere_qv[0, :, :, :]
+            _val_init_atmosphere_u[:, :, :] = init_atmosphere_u[0, :, :, 1:]
+            _val_init_atmosphere_v[:, :, :] = init_atmosphere_v[0, :, 1:, :]
+            _val_init_atmosphere_w[:, :, :] = init_atmosphere_w[0, :, :, :]
+            _val_init_soil_t[:, :, :] = init_soil_t[0, :, :, :]
+            for k in range(0,_val_init_soil_m.shape[0]):
+                # adjust soil moisture according soil_moisture_adjust field (if exists)
+                _val_init_soil_m[k, :, :] = init_soil_m[0, k, :, :] * sma[:, :]
+        finally:
+            # close interpolated file
+            infile.close()
+        #
+        if not nested_domain:
+            # cycle over all included time steps
+            for ts in range(0, len(interp_files)):
+                print("Processing file",interp_files[ts])
+                # open corresponding interpolated file
+                infile = netCDF4.Dataset(interp_files[ts], "r", format="NETCDF4")
+                try:
+                    # open variables in the input file
+                    init_atmosphere_pt = infile.variables['init_atmosphere_pt']
+                    init_atmosphere_qv = infile.variables['init_atmosphere_qv']
+                    init_atmosphere_u = infile.variables['init_atmosphere_u']
+                    init_atmosphere_v = infile.variables['init_atmosphere_v']
+                    init_atmosphere_w = infile.variables['init_atmosphere_w']
+                    surface_forcing_surface_pressure = infile.variables['surface_forcing_surface_pressure']
+
+                    ##################################
+                    # write values for time dependent values
+                    # surface pressure
+                    _val_surface_forcing_surface_pressure[ts] = np.average(surface_forcing_surface_pressure[:,:,:], axis = (1,2))[0]
+                    # boundary conditions
+                    _val_ls_forcing_pt_left[ts, :, :] = init_atmosphere_pt[0, :, :, 0]
+                    _val_ls_forcing_pt_right[ts, :, :] = init_atmosphere_pt[0, :, :, dimensions['xdim'] - 1]
+                    _val_ls_forcing_pt_south[ts, :, :] = init_atmosphere_pt[0, :, 0, :]
+                    _val_ls_forcing_pt_north[ts, :, :] = init_atmosphere_pt[0, :, dimensions['ydim'] - 1, :]
+                    _val_ls_forcing_pt_top[ts, :, :] = init_atmosphere_pt[0, dimensions['zdim'] - 1, :, :]
+
+                    _val_ls_forcing_qv_left[ts, :, :] = init_atmosphere_qv[0, :, :, 0]
+                    _val_ls_forcing_qv_right[ts, :, :] = init_atmosphere_qv[0, :, :, dimensions['xdim'] - 1]
+                    _val_ls_forcing_qv_south[ts, :, :] = init_atmosphere_qv[0, :, 0, :]
+                    _val_ls_forcing_qv_north[ts, :, :] = init_atmosphere_qv[0, :, dimensions['ydim'] - 1, :]
+                    _val_ls_forcing_qv_top[ts, :, :] = init_atmosphere_qv[0, dimensions['zdim'] - 1, :, :]
+
+                    # Perform mass balancing
+                    uxleft = init_atmosphere_u[0, :, :, 0]
+                    uxright = init_atmosphere_u[0, :, :, dimensions['xdim'] - 1]
+                    vysouth = init_atmosphere_v[0, :, 0, :]
+                    vynorth = init_atmosphere_v[0, :, dimensions['ydim'] - 1, :]
+                    wztop = init_atmosphere_w[0, dimensions['zwdim'] - 1, :, :]
+                    mass_disbalance = ((uxleft * areas_xb).sum()
+                        - (uxright * areas_xb).sum()
+                        + (vysouth * areas_yb).sum()
+                        - (vynorth * areas_yb).sum()
+                        - (wztop * areas_zb).sum())
+                    mass_corr_v = mass_disbalance / area_boundaries
+                    print('Mass disbalance: {0:8g} m3/s (avg = {1:8g} m/s)'.format(
+                        mass_disbalance, mass_corr_v))
+                    uxleft -= mass_corr_v
+                    uxright += mass_corr_v
+                    vysouth -= mass_corr_v
+                    vynorth += mass_corr_v
+                    wztop += mass_corr_v
+
+                    # Verify mass balance (optional)
+                    #mass_disbalance = ((uxleft * areas_xb).sum()
+                    #    - (uxright * areas_xb).sum()
+                    #    + (vysouth * areas_yb).sum()
+                    #    - (vynorth * areas_yb).sum()
+                    #    - (wztop * areas_zb).sum())
+                    #mass_corr_v = mass_disbalance / area_boundaries
+                    #print('Mass balanced:   {0:8g} m3/s (avg = {1:8g} m/s)'.format(
+                    #    mass_disbalance, mass_corr_v))
+
+                    _val_ls_forcing_u_left[ts, :, :] = uxleft
+                    _val_ls_forcing_u_right[ts, :, :] = uxright
+                    _val_ls_forcing_u_south[ts, :, :] = init_atmosphere_u[0, :, 0, 1:]
+                    _val_ls_forcing_u_north[ts, :, :] = init_atmosphere_u[0, :, dimensions['ydim'] - 1, 1:]
+                    _val_ls_forcing_u_top[ts, :, :] = init_atmosphere_u[0, dimensions['zdim'] - 1, :, 1:]
+
+                    _val_ls_forcing_v_left[ts, :, :] = init_atmosphere_v[0, :, 1:, 0]
+                    _val_ls_forcing_v_right[ts, :, :] = init_atmosphere_v[0, :, 1:, dimensions['xdim'] - 1]
+                    _val_ls_forcing_v_south[ts, :, :] = vysouth
+                    _val_ls_forcing_v_north[ts, :, :] = vynorth
+                    _val_ls_forcing_v_top[ts, :, :] = init_atmosphere_v[0, dimensions['zdim'] - 1, 1:, :]
+
+                    _val_ls_forcing_w_left[ts, :, :] = init_atmosphere_w[0, :, :, 0]
+                    _val_ls_forcing_w_right[ts, :, :] = init_atmosphere_w[0, :, :, dimensions['xdim'] - 1]
+                    _val_ls_forcing_w_south[ts, :, :] = init_atmosphere_w[0, :, 0, :]
+                    _val_ls_forcing_w_north[ts, :, :] = init_atmosphere_w[0, :, dimensions['ydim'] - 1, :]
+                    _val_ls_forcing_w_top[ts, :, :] = wztop
+
+                finally:
+                    # close interpolated file
+                    infile.close()
+
+                # write geostrophic wind
+                print('Open wrf file '+wrf_files[ts])
+                nc_wrf = netCDF4.Dataset(wrf_files[ts], 'r')
+                try:
+                    ug, vg = palm_wrf_gw(nc_wrf, lon_center, lat_center, z_levels)
+                    _val_ls_forcing_ug[ts, :] = ug
+                    _val_ls_forcing_vg[ts, :] = vg
+                finally:
+                    nc_wrf.close()
+
+            # Write chemical boundary conds
+            if camx_interp_fname:
+                f_camx = netCDF4.Dataset(camx_interp_fname)
+                try:
+                    for vname, vval in f_camx.variables.items():
+                        # PALM doesn't support 3D LOD=2 init for chem yet, we have to average the field
+                        var = outfile.createVariable('init_atmosphere_'+vname,
+                                'f4', ('z',), fill_value=fillvalue_float)
+                        var.units = vval.units
+                        var.lod = 1
+                        var[:] = vval[0,:,:,:].mean(axis=(1,2))
+
+                        var = outfile.createVariable('ls_forcing_left_'+vname,
+                                'f4', ('time','z','y'), fill_value=fillvalue_float)
+                        var.units = vval.units
+                        var[:] = vval[:,:,:,0]
+
+                        var = outfile.createVariable('ls_forcing_right_'+vname,
+                                'f4', ('time','z','y'), fill_value=fillvalue_float)
+                        var.units = vval.units
+                        var[:] = vval[:,:,:,-1]
+
+                        var = outfile.createVariable('ls_forcing_south_'+vname,
+                                'f4', ('time','z','x'), fill_value=fillvalue_float)
+                        var.units = vval.units
+                        var[:] = vval[:,:,0,:]
+
+                        var = outfile.createVariable('ls_forcing_north_'+vname,
+                                'f4', ('time','z','x'), fill_value=fillvalue_float)
+                        var.units = vval.units
+                        var[:] = vval[:,:,-1,:]
+
+                        var = outfile.createVariable('ls_forcing_top_'+vname,
+                                'f4', ('time','y','x'), fill_value=fillvalue_float)
+                        var.units = vval.units
+                        var[:] = vval[:,-1,:,:]
+                finally:
+                    f_camx.close()
+
+        if len(rad_times_proc) > 0:
+            # process radiation inputs
+            # radiation time dimension and variable
+            outfile.createDimension('time_rad', len(rad_times_proc))
+            _val_times =  outfile.createVariable('time_rad',"f4", ("time_rad"))
+            _val_times[:] = rad_times_proc[:]
+            # radiation variables
+            var = outfile.createVariable('rad_sw_in', "f4", ("time_rad"), fill_value=fillvalue_float)
+            var.setncattr('lod', 1)
+            var.units = 'W/m2'
+            var[:] = rad_values_proc[0][:]
+            var = outfile.createVariable('rad_lw_in', "f4", ("time_rad"), fill_value=fillvalue_float)
+            var.setncattr('lod', 1)
+            var.units = 'W/m2'
+            var[:] = rad_values_proc[1][:]
+            var = outfile.createVariable('rad_sw_in_dif', "f4", ("time_rad"), fill_value=fillvalue_float)
+            var.setncattr('lod', 1)
+            var.units = 'W/m2'
+            var[:] = rad_values_proc[2][:]
+
+    finally:
+        outfile.close()
+    
+
+
+
+
Index: /palm/trunk/UTIL/WRF_interface/dynamic/palm_wrf_utils.py
===================================================================
--- /palm/trunk/UTIL/WRF_interface/dynamic/palm_wrf_utils.py	(revision 4766)
+++ /palm/trunk/UTIL/WRF_interface/dynamic/palm_wrf_utils.py	(revision 4766)
@@ -0,0 +1,643 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+#------------------------------------------------------------------------------#
+#
+# Scripts for processing of WRF and CAMx files to PALM dynamic driver
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+#
+# Copyright 2018-2020 Institute of Computer Science
+#                     of the Czech Academy of Sciences, Prague
+# Authors: Krystof Eben, Jaroslav Resler, Pavel Krc
+#
+#------------------------------------------------------------------------------#
+
+'''WRF (+CAMx) utility module for PALM dynamic driver generator'''
+
+import os
+import math
+import numpy as np
+import pyproj
+import scipy.ndimage as ndimage
+import netCDF4
+
+import metpy.calc as mpcalc
+from metpy.interpolate import interpolate_1d, log_interpolate_1d
+from metpy.units import units
+
+from palm_dynamic_config import *
+
+# Constants directly equivalent to WRF code
+radius = 6370000.0
+g = 9.81 #m/s2
+rd = 287. #dry air gas constant (J/kg/K)
+rd_cp = 2./7. #from WRF v4 technote (R_d / c_p)
+wrf_base_temp = 300. #NOT wrfout T00
+
+_ax = np.newaxis
+
+class WRFCoordTransform(object):
+    'Coordinate transformer for WRFOUT files'
+
+    def __init__(self, ncf):
+        attr = lambda a: getattr(ncf, a)
+
+        # Define grids
+
+        # see http://www.pkrc.net/wrf-lambert.html
+        #latlon_wgs84 = pyproj.Proj(proj='latlong',
+        #    ellps='WGS84', datum='WGS84',
+        #    no_defs=True) #don't use - WRF datum misuse
+
+        latlon_sphere = pyproj.Proj(proj='latlong',
+            a=radius, b=radius,
+            towgs84='0,0,0', no_defs=True)
+
+        lambert_grid = pyproj.Proj(proj='lcc',
+            lat_1=attr('TRUELAT1'),
+            lat_2=attr('TRUELAT2'),
+            lat_0=attr('MOAD_CEN_LAT'),
+            lon_0=attr('STAND_LON'),
+            a=radius, b=radius,
+            towgs84='0,0,0', no_defs=True)
+
+        # resoltion in m
+        self.dx = dx = attr('DX')
+        self.dy = dy = attr('DY')
+
+        # number of mass grid points
+        self.nx = nx = attr('WEST-EAST_GRID_DIMENSION') - 1
+        self.ny = ny = attr('SOUTH-NORTH_GRID_DIMENSION') - 1
+
+        # distance between centers of mass grid points at edges
+        extent_x = (nx - 1) * dx
+        extent_y = (ny - 1) * dy
+
+        # grid center in lambert
+        center_x, center_y = pyproj.transform(latlon_sphere, lambert_grid,
+            attr('CEN_LON'), attr('CEN_LAT'))
+
+        # grid origin coordinates in lambert
+        i0_x = center_x - extent_x*.5
+        j0_y = center_y - extent_y*.5
+
+        # Define fast transformation methods
+
+        def latlon_to_ji(lat, lon):
+            x, y = pyproj.transform(latlon_sphere, lambert_grid,
+                    lon, lat)
+            return (y-j0_y)/dy, (x-i0_x)/dx
+        self.latlon_to_ji = latlon_to_ji
+
+        def ji_to_latlon(j, i):
+            lon, lat = pyproj.transform(lambert_grid, latlon_sphere,
+                i*dx+i0_x, j*dy+j0_y)
+            return lat, lon
+        self.ji_to_latlon = ji_to_latlon
+
+    def verify(self, ncf):
+        lat = ncf.variables['XLAT'][0]
+        lon = ncf.variables['XLONG'][0]
+        j, i = np.mgrid[0:self.ny, 0:self.nx]
+
+        jj, ii = self.latlon_to_ji(lat, lon)
+        d = np.hypot(jj-j, ii-i)
+        print('error for ll->ji: max {0} m, avg {1} m.'.format(d.max(), d.mean()))
+
+        llat, llon = self.ji_to_latlon(j, i)
+        d = np.hypot(llat - lat, llon - lon)
+        print('error for ji->ll: max {0} deg, avg {1} deg.'.format(d.max(), d.mean()))
+
+        lat = ncf.variables['XLAT_U'][0]
+        lon = ncf.variables['XLONG_U'][0]
+        j, i = np.mgrid[0:self.ny, 0:self.nx+1]
+        jj, ii = self.latlon_to_ji(lat, lon)
+        ii = ii + .5
+        d = np.hypot(jj-j, ii-i)
+        print('error for U-staggered ll->ji: max {0} m, avg {1} m.'.format(d.max(), d.mean()))
+
+class CAMxCoordTransform(object):
+    'Coordinate transformer for CAMx files running from WRF'
+
+    def __init__(self, ncf):
+        attr = lambda a: getattr(ncf, a)
+
+        # Define grids
+
+        latlon_sphere = pyproj.Proj(proj='latlong',
+            a=radius, b=radius,
+            towgs84='0,0,0', no_defs=True)
+
+        lambert_grid = pyproj.Proj(proj='lcc',
+            lat_1=attr('P_ALP'),
+            lat_2=attr('P_BET'),
+            lat_0=attr('YCENT'),
+            lon_0=attr('P_GAM'),
+            a=radius, b=radius,
+            towgs84='0,0,0', no_defs=True)
+
+        # resoltion in m
+        self.dx = dx = attr('XCELL')
+        self.dy = dy = attr('YCELL')
+
+        # number of mass grid points
+        self.nx = nx = attr('NCOLS')
+        self.ny = ny = attr('NROWS')
+
+        # grid origin coordinates in lambert
+        i0_x = attr('XORIG')
+        j0_y = attr('YORIG')
+
+        # Define fast transformation methods
+
+        def latlon_to_ji(lat, lon):
+            x, y = pyproj.transform(latlon_sphere, lambert_grid,
+                    lon, lat)
+            return (y-j0_y)/dy, (x-i0_x)/dx
+        self.latlon_to_ji = latlon_to_ji
+
+        def ji_to_latlon(j, i):
+            lon, lat = pyproj.transform(lambert_grid, latlon_sphere,
+                i*dx+i0_x, j*dy+j0_y)
+            return lat, lon
+        self.ji_to_latlon = ji_to_latlon
+
+    def verify(self, ncf):
+        lat = ncf.variables['latitude'][:]
+        lon = ncf.variables['longitude'][:]
+        j, i = np.mgrid[0:self.ny, 0:self.nx]
+
+        jj, ii = self.latlon_to_ji(lat, lon)
+        d = np.hypot(jj-j, ii-i)
+        print('error for ll->ji: max {0} m, avg {1} m.'.format(d.max(), d.mean()))
+
+        llat, llon = self.ji_to_latlon(j, i)
+        d = np.hypot(llat - lat, llon - lon)
+        print('error for ji->ll: max {0} deg, avg {1} deg.'.format(d.max(), d.mean()))
+
+class BilinearRegridder(object):
+    '''Bilinear regridder for multidimensional data.
+
+    By standard, the last two dimensions are always Y,X in that order.
+    '''
+    def __init__(self, projected_x, projected_y, preloaded=False):
+        projected_x = np.asanyarray(projected_x)
+        projected_y = np.asanyarray(projected_y)
+        self.shape = projected_x.shape
+        self.rank = len(self.shape)
+        assert self.shape == projected_y.shape
+
+        y0 = np.floor(projected_y)
+        yd = projected_y - y0
+        ydd = 1. - yd
+        self.y0 = y0.astype('i8')
+
+        x0 = np.floor(projected_x)
+        xd = projected_x - x0
+        xdd = 1. - xd
+        self.x0 = x0.astype('i8')
+
+        if preloaded:
+            # Prepare slices for preloading from NetCDF files (for cases where
+            # the range of loaded Y, X coordinates is much less than total
+            # size. The regrid method then expects preloaded data.
+
+            ybase = self.y0.min()
+            self.ys = slice(ybase, self.y0.max()+2)
+            self.y0 -= ybase
+
+            xbase = self.x0.min()
+            self.xs = slice(xbase, self.x0.max()+2)
+            self.x0 -= xbase
+
+        self.y1 = self.y0 + 1
+        self.x1 = self.x0 + 1
+
+        self.weights = np.array([
+            ydd * xdd, #wy0x0
+            ydd * xd , #wy0x1
+            yd  * xdd, #wy1x0
+            yd  * xd , #wy1x1
+            ])
+
+    def regrid(self, data):
+        # data may contain additional dimensions (before Y,X)
+        dshape = data.shape[:-2]
+        drank = len(dshape)
+
+        # Prepare array for selected data
+        sel_shape = (4,) + dshape + self.shape
+        selection = np.empty(sel_shape, dtype=data.dtype)
+
+        selection[0, ...] = data[..., self.y0, self.x0]
+        selection[1, ...] = data[..., self.y0, self.x1]
+        selection[2, ...] = data[..., self.y1, self.x0]
+        selection[3, ...] = data[..., self.y1, self.x1]
+
+        # Slice weights to match the extra dimensions
+        wslice = ((slice(None),) +      #weights
+            (np.newaxis,) * drank +     #data minus Y,X
+            (slice(None),) * self.rank) #regridded shape
+
+        w = selection * self.weights[wslice]
+        return w.sum(axis=0)
+
+def print_dstat(desc, delta):
+    print('Delta stats for {0} ({1:8g} ~ {2:8g}): bias = {3:8g}, MAE = {4:8g}, RMSE = {5:8g}'.format(
+        desc, delta.min(), delta.max(), delta.mean(), np.abs(delta).mean(), np.sqrt((delta**2).mean())))
+
+def barom_pres(p0, gp, gp0, t0):
+    barom = 1. / (rd * t0)
+    return p0 * np.exp((gp0-gp)*barom)
+
+def barom_gp(gp0, p, p0, t0):
+    baromi = rd * t0
+    return gp0 - np.log(p/p0) * baromi
+
+def calc_ph_hybrid(f, mu):
+    pht = f.variables['P_TOP'][0]
+    c3f = f.variables['C3F'][0]
+    c4f = f.variables['C4F'][0]
+    c3h = f.variables['C3H'][0]
+    c4h = f.variables['C4H'][0]
+    return (c3f[:,_ax,_ax]*mu[_ax,:,:] + (c4f[:,_ax,_ax] + pht),
+            c3h[:,_ax,_ax]*mu[_ax,:,:] + (c4h[:,_ax,_ax] + pht))
+
+def calc_ph_sigma(f, mu):
+    pht = f.variables['P_TOP'][0]
+    eta_f = f.variables['ZNW'][0]
+    eta_h = f.variables['ZNU'][0]
+    return (eta_f[:,_ax,_ax]*mu[_ax,:,:] + pht,
+            eta_h[:,_ax,_ax]*mu[_ax,:,:] + pht)
+
+def wrf_t(f):
+    p = f.variables['P'][0,:,:,:] + f.variables['PB'][0,:,:,:]
+    return (f.variables['T'][0,:,:,:] + wrf_base_temp) * np.power(0.00001*p, rd_cp)
+
+def calc_gp(f, ph):
+    terr = f.variables['HGT'][0,:,:]
+    gp0 = terr * g
+    gp = [gp0]
+    t = wrf_t(f)
+    for lev in range(1, ph.shape[0]):
+        gp.append(barom_gp(gp[-1], ph[lev,:,:], ph[lev-1,:,:], t[lev-1,:,:]))
+    return np.array(gp)
+
+def palm_wrf_vertical_interp(infile, outfile, wrffile, z_levels, z_levels_stag,
+        z_soil_levels, origin_z, terrain, wrf_hybrid_levs, vinterp_terrain_smoothing):
+
+    zdim = len(z_levels)
+    zwdim = len(z_levels_stag)
+    zsoildim = len(z_soil_levels)
+    dimnames = ['z', 'zw', 'zsoil']   # dimnames of palm vertical dims
+    dimensions = [zdim , zwdim, zsoildim]
+
+    print("infile: " , infile)
+    print("outfile: ", outfile)
+
+    try:
+        os.remove(outfile)
+        os.remove(infile+'_vinterp.log')
+    except:
+        pass
+
+    nc_infile = netCDF4.Dataset(infile, 'r')
+    nc_wrf = netCDF4.Dataset(wrffile, 'r')
+    nc_outfile = netCDF4.Dataset(outfile, "w", format="NETCDF4")
+    nc_outfile.createDimension('Time', None)
+    for dimname in ['west_east', 'south_north', 'soil_layers_stag']:
+        nc_outfile.createDimension(dimname, len(nc_infile.dimensions[dimname]))
+    for i in range (len(dimnames)):
+        nc_outfile.createDimension(dimnames[i], dimensions[i])
+
+    # Use hybrid ETA levels in WRF and stretch them so that the WRF terrain
+    # matches either PALM terrain or flat terrain at requested height
+    gpf = nc_infile.variables['PH'][0,:,:,:] + nc_infile.variables['PHB'][0,:,:,:]
+    wrfterr = gpf[0]*(1./g)
+
+    if vinterp_terrain_smoothing is None:
+        target_terrain = terrain
+    else:
+        print('Smoothing PALM terrain for the purpose of dynamic driver with sigma={0} grid points.'.format(
+            vinterp_terrain_smoothing))
+        target_terrain = ndimage.gaussian_filter(terrain, sigma=vinterp_terrain_smoothing, order=0)
+    print('Morphing WRF terrain ({0} ~ {1}) to PALM terrain ({2} ~ {3})'.format(
+        wrfterr.min(), wrfterr.max(), target_terrain.min(), target_terrain.max()))
+    print_dstat('terrain shift', wrfterr - target_terrain[:,:])
+
+    # Load original dry air column pressure
+    mu = nc_infile.variables['MUB'][0,:,:] + nc_infile.variables['MU'][0,:,:]
+    pht = nc_wrf.variables['P_TOP'][0]
+
+    # Shift column pressure so that it matches PALM terrain
+    t = wrf_t(nc_infile)
+    mu2 = barom_pres(mu+pht, target_terrain*g, gpf[0,:,:], t[0,:,:])-pht
+
+    # Calculate original and shifted 3D dry air pressure
+    if wrf_hybrid_levs:
+        phf, phh = calc_ph_hybrid(nc_wrf, mu)
+        phf2, phh2 = calc_ph_hybrid(nc_wrf, mu2)
+    else:
+        phf, phh = calc_ph_sigma(nc_wrf, mu)
+        phf2, phh2 = calc_ph_sigma(nc_wrf, mu2)
+
+    # Shift 3D geopotential according to delta dry air pressure
+    tf = np.concatenate((t, t[-1:,:,:]), axis=0) # repeat highest layer
+    gpf2 = barom_gp(gpf, phf2, phf, tf)
+    # For half-levs, originate from gp full levs rather than less accurate gp halving
+    gph2 = barom_gp(gpf[:-1,:,:], phh2, phf[:-1,:,:], t)
+    zf = gpf2 * (1./g) - origin_z
+    zh = gph2 * (1./g) - origin_z
+
+    # Report
+    gpdelta = gpf2 - gpf
+    print('GP deltas by level:')
+    for k in range(gpf.shape[0]):
+        print_dstat(k, gpdelta[k])
+
+    # Because we require levels below the lowest level from WRF, we will always
+    # add one layer at zero level with repeated values from the lowest level.
+    # WRF-python had some special treatment for theta in this case.
+    height = np.zeros((zh.shape[0]+1,) + zh.shape[1:], dtype=zh.dtype)
+    height[0,:,:] = -999. #always below terrain
+    height[1:,:,:] = zh
+    heightw = np.zeros((zf.shape[0]+1,) + zf.shape[1:], dtype=zf.dtype)
+    heightw[0,:,:] = -999. #always below terrain
+    heightw[1:,:,:] = zf
+
+    # ======================== SPECIFIC HUMIDITY ==============================
+    qv_raw = nc_infile.variables['SPECHUM'][0]
+    qv_raw = np.r_[qv_raw[0:1], qv_raw]
+    
+    # Vertical interpolation to grid height levels (specified in km!)
+    # Levels start at 50m (below that the interpolation looks very sketchy)
+    init_atmosphere_qv = interpolate_1d(z_levels, height, qv_raw)
+    vdata = nc_outfile.createVariable('init_atmosphere_qv', "f4", ("Time", "z","south_north","west_east"))
+    vdata[0,:,:,:] = init_atmosphere_qv
+    
+    # ======================= POTENTIAL TEMPERATURE ==========================
+    pt_raw = nc_infile.variables['T'][0] + 300.   # from perturbation pt to standard
+    pt_raw = np.r_[pt_raw[0:1], pt_raw]
+    init_atmosphere_pt = interpolate_1d(z_levels, height, pt_raw)
+    
+    #plt.figure(); plt.contourf(pt[0]) ; plt.colorbar() ; plt.show()
+    vdata = nc_outfile.createVariable('init_atmosphere_pt', "f4", ("Time", "z","south_north","west_east"))
+    vdata[0,:,:,:] = init_atmosphere_pt
+    
+    # ======================= Wind ==========================================
+    u_raw = nc_infile.variables['U'][0]
+    u_raw = np.r_[u_raw[0:1], u_raw]
+    init_atmosphere_u = interpolate_1d(z_levels, height, u_raw)
+    
+    vdata = nc_outfile.createVariable('init_atmosphere_u', "f4", ("Time", "z","south_north","west_east"))
+    vdata[0,:,:,:] = init_atmosphere_u
+
+    v_raw = nc_infile.variables['V'][0]
+    v_raw = np.r_[v_raw[0:1], v_raw]
+    init_atmosphere_v = interpolate_1d(z_levels, height, v_raw)
+    
+    vdata = nc_outfile.createVariable('init_atmosphere_v', "f4", ("Time", "z","south_north","west_east"))
+    #vdata.coordinates = "XLONG_V XLAT_V XTIME"
+    vdata[0,:,:,:] = init_atmosphere_v
+    
+    w_raw = nc_infile.variables['W'][0]
+    w_raw = np.r_[w_raw[0:1], w_raw]
+    init_atmosphere_w = interpolate_1d(z_levels_stag, heightw, w_raw)
+    
+    vdata = nc_outfile.createVariable('init_atmosphere_w', "f4", ("Time", "zw","south_north","west_east"))
+    #vdata.coordinates = "XLONG XLAT XTIME"
+    vdata[0,:,:,:] = init_atmosphere_w
+
+    # ===================== SURFACE PRESSURE ==================================
+    surface_forcing_surface_pressure = nc_infile.variables['PSFC']
+    vdata = nc_outfile.createVariable('surface_forcing_surface_pressure', "f4", ("Time", "south_north","west_east"))
+    vdata[0,:,:] = surface_forcing_surface_pressure[0,:,:]
+    
+    # ======================== SOIL VARIABLES (without vertical interpolation) =============
+    # soil temperature
+    init_soil_t = nc_infile.variables['TSLB']
+    vdata = nc_outfile.createVariable('init_soil_t', "f4", ("Time", "zsoil","south_north","west_east"))
+    vdata[0,:,:,:] = init_soil_t[0,:,:,:]
+
+    # soil moisture
+    init_soil_m = nc_infile.variables['SMOIS']
+    vdata = nc_outfile.createVariable('init_soil_m', "f4", ("Time","zsoil","south_north","west_east"))
+    vdata[0,:,:,:] = init_soil_m[0,:,:,:]
+
+    # zsoil
+    zsoil = nc_wrf.variables['ZS']    #ZS:description = "DEPTHS OF CENTERS OF SOIL LAYERS" ;
+    vdata = nc_outfile.createVariable('zsoil', "f4", ("zsoil"))
+    vdata[:] = zsoil[0,:]
+
+    # coordinates z, zw
+    vdata = nc_outfile.createVariable('z', "f4", ("z"))
+    vdata[:] = list(z_levels)
+
+    vdata = nc_outfile.createVariable('zw', "f4", ("zw"))
+    vdata[:] = list (z_levels_stag) 
+
+    # zsoil is taken from wrf - not need to define it
+
+    nc_infile.close()
+    nc_wrf.close()
+    nc_outfile.close()
+
+def palm_wrf_gw(f, lon, lat, levels, tidx=0):
+    '''Calculate geostrophic wind from WRF using metpy'''
+
+    hgts, ug, vg = calcgw_wrf(f, lat, lon, levels, tidx)
+
+    # extrapolate at the bottom
+    hgts = np.r_[np.array([0.]), hgts]
+    ug = np.r_[ug[0], ug]
+    vg = np.r_[vg[0], vg]
+
+    return minterp(levels, hgts, ug, vg)
+
+
+def minterp(interp_heights, data_heights, u, v):
+    '''Interpolate wind using power law for agl levels'''
+
+    pdata = data_heights ** gw_alpha
+    pinterp = interp_heights ** gw_alpha
+    hindex = np.searchsorted(data_heights, interp_heights, side='right')
+    lindex = hindex - 1
+    assert lindex[0] >= 0
+    assert hindex[-1] < len(data_heights)
+    lbound = pdata[lindex]
+    hcoef = (pinterp - lbound) / (pdata[hindex] - lbound)
+    #print(data_heights)
+    #print(lindex)
+    #print(hcoef)
+    lcoef = 1. - hcoef
+    iu = u[lindex] * lcoef + u[hindex] * hcoef
+    iv = v[lindex] * lcoef + v[hindex] * hcoef
+    return iu, iv
+
+def get_wrf_dims(f, lat, lon, xlat, xlong):
+    '''A crude method, yet satisfactory for approximate WRF surroundings'''
+
+    sqdist = (xlat - lat)**2 + (xlong - lon)**2
+    coords = np.unravel_index(sqdist.argmin(), sqdist.shape)
+
+    xmargin = int(math.ceil(gw_wrf_margin_km * 1000 / f.DX)) #py2 ceil produces float
+    ymargin = int(math.ceil(gw_wrf_margin_km * 1000 / f.DY))
+    y0, y1 = coords[0] - ymargin, coords[0] + ymargin
+    x0, x1 = coords[1] - xmargin, coords[1] + xmargin
+    assert 0 <= y0 < y1 < sqdist.shape[0], "Point {0} + surroundings not inside domain".format(coords[0])
+    assert 0 <= x0 < x1 < sqdist.shape[1], "Point {0} + surroundings not inside domain".format(coords[1])
+
+    return coords, (slice(y0, y1+1), slice(x0, x1+1)), (ymargin, xmargin)
+
+def calcgw_wrf(f, lat, lon, levels, tidx=0):
+    # MFDataset removes the time dimension from XLAT, XLONG
+    xlat = f.variables['XLAT']
+    xlslice = (0,) * (len(xlat.shape)-2) + (slice(None), slice(None))
+    xlat = xlat[xlslice]
+    xlong = f.variables['XLONG'][xlslice]
+
+    (iy, ix), area, (iby, ibx) = get_wrf_dims(f, lat, lon, xlat, xlong)
+    areat = (tidx,) + area
+    areatz = (tidx, slice(None)) + area
+    #print('wrf coords', lat, lon, xlat[iy,ix], xlong[iy,ix])
+    #print(xlat[area][iby,ibx], xlong[area][iby,ibx], areat)
+
+    # load area
+    hgt = (f.variables['PH'][areatz] + f.variables['PHB'][areatz]) / 9.81
+    hgtu = (hgt[:-1] + hgt[1:]) * .5
+    pres = f.variables['P'][areatz] + f.variables['PB'][areatz]
+    terrain = f.variables['HGT'][areat]
+
+    # find suitable pressure levels
+    yminpres, xminpres = np.unravel_index(pres[0].argmin(), pres[0].shape)
+    pres1 = pres[0, yminpres, xminpres] - 1.
+
+    aglpt = hgtu[:,iby,ibx] - terrain[iby,ibx]
+    pres0 = pres[np.searchsorted(aglpt, levels[-1]), iby, ibx]
+    plevels = np.arange(pres1, min(pres0, pres1)-1, -1000.)
+
+    # interpolate wrf into pressure levels
+    phgt = log_interpolate_1d(plevels, pres, hgtu, axis=0)
+
+    # Set up some constants based on our projection, including the Coriolis parameter and
+    # grid spacing, converting lon/lat spacing to Cartesian
+    coriol = mpcalc.coriolis_parameter(np.deg2rad(xlat[area])).to('1/s')
+
+    # lat_lon_grid_deltas doesn't work under py2, but for WRF grid it is still
+    # not very accurate, better use direct values.
+    #dx, dy = mpcalc.lat_lon_grid_deltas(xlong[area], xlat[area])
+    dx = f.DX * units.m
+    dy = f.DY * units.m
+
+    # Smooth height data. Sigma=1.5 for gfs 0.5deg
+    res_km = f.DX / 1000.
+
+    ug = np.zeros(plevels.shape, 'f8')
+    vg = np.zeros(plevels.shape, 'f8')
+    for i in range(len(plevels)):
+        sh = ndimage.gaussian_filter(phgt[i,:,:], sigma=1.5*50/res_km, order=0)
+        ugl, vgl = mpcalc.geostrophic_wind(sh * units.m, coriol, dx, dy)
+        ug[i] = ugl[iby, ibx].magnitude
+        vg[i] = vgl[iby, ibx].magnitude
+
+    return phgt[:,iby,ibx], ug, vg
+
+# The following two functions calculate GW from GFS files, although this
+# function is currently not implemented in PALM dynamic driver generation
+# script
+
+def calcgw_gfs(v, lat, lon):
+    height, lats, lons = v.data(lat1=lat-gw_gfs_margin_deg ,lat2=lat+gw_gfs_margin_deg,
+            lon1=lon-gw_gfs_margin_deg, lon2=lon+gw_gfs_margin_deg)
+    i = np.searchsorted(lats[:,0], lat)
+    if abs(lats[i+1,0] - lat) < abs(lats[i,0] - lat):
+        i = i+1
+    j = np.searchsorted(lons[0,:], lon)
+    if abs(lons[0,i+1] - lon) < abs(lons[0,i] - lon):
+        j = j+1
+    #print('level', v.level, 'height', height[i,j], lats[i,j], lons[i,j])
+
+    # Set up some constants based on our projection, including the Coriolis parameter and
+    # grid spacing, converting lon/lat spacing to Cartesian
+    f = mpcalc.coriolis_parameter(np.deg2rad(lats)).to('1/s')
+    dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats)
+    res_km = (dx[i,j]+dy[i,j]).magnitude / 2000.
+
+    # Smooth height data. Sigma=1.5 for gfs 0.5deg
+    height = ndimage.gaussian_filter(height, sigma=1.5*50/res_km, order=0)
+
+    # In MetPy 0.5, geostrophic_wind() assumes the order of the dimensions is (X, Y),
+    # so we need to transpose from the input data, which are ordered lat (y), lon (x).
+    # Once we get the components,transpose again so they match our original data.
+    geo_wind_u, geo_wind_v = mpcalc.geostrophic_wind(height * units.m, f, dx, dy)
+
+    return height[i,j], geo_wind_u[i,j], geo_wind_v[i,j]
+
+def combinegw_gfs(grbs, levels, lat, lon):
+    heights = []
+    us = []
+    vs = []
+    for grb in grbs:
+        h, u, v = calcgw_gfs(grb, lat, lon)
+        heights.append(h)
+        us.append(u.magnitude)
+        vs.append(v.magnitude)
+    heights = np.array(heights)
+    us = np.array(us)
+    vs = np.array(vs)
+
+    ug, vg = minterp(np.asanyarray(levels), heights[::-1], us[::-1], vs[::-1])
+    return ug, vg
+
+if __name__ == '__main__':
+    import argparse
+
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument('-w', '--wrfout', help='verify wrfout file')
+    parser.add_argument('-c', '--camx', help='verify camx file')
+    args = parser.parse_args()
+
+    if args.wrfout:
+        f = netCDF4.Dataset(args.wrfout)
+
+        print('Verifying coord transform:')
+        t = WRFCoordTransform(f)
+        t.verify(f)
+
+        print('\nVerifying vertical levels:')
+        mu = f.variables['MUB'][0,:,:] + f.variables['MU'][0,:,:]
+        gp = f.variables['PH'][0,:,:,:] + f.variables['PHB'][0,:,:,:]
+
+        print('\nUsing sigma:')
+        phf, phh = calc_ph_sigma(f, mu)
+        gp_calc = calc_gp(f, phf)
+        delta = gp_calc - gp
+        for lev in range(delta.shape[0]):
+            print_dstat(lev, delta[lev])
+
+        print('\nUsing hybrid:')
+        phf, phh = calc_ph_hybrid(f, mu)
+        gp_calc = calc_gp(f, phf)
+        delta = gp_calc - gp
+        for lev in range(delta.shape[0]):
+            print_dstat(lev, delta[lev])
+
+        f.close()
+
+    if args.camx:
+        f = netCDF4.Dataset(args.camx)
+        t = CAMxCoordTransform(f)
+        t.verify(f)
+        f.close()
