def create_librad_netcdf_cldfile(fname, lwc, reff, hhl, dx, dy, reff_min=2.51, reff_max=55):
    """
    Creates a netcdf cloud file for libRadtran calculations
    Provide:
    lwc and reff arrays with dimensions (Nx, Ny, Nz)
    hhl with dimension(Nz+1)

    dx, dy, hhl in units of [km]
    lwc in [g/kg]
    reff in [micrometer]
    """
    import xarray as xr
    import os
    import numpy as np
    #pylint: disable = invalid-name
    """ create lwc file """
    if (lwc is None) or (reff is None):
      print('lwc or reff is None, not writing a cld file...')
      return

    print('shape of input:', np.shape(lwc,), ':', np.shape(reff), ':', np.shape(hhl))
    print('assembling data for:  ', fname, ' :: min {} max {}'.format(np.min(lwc), np.max(lwc)), ' :: reff {} {}'.format(np.min(reff), np.max(reff)))
    if lwc.shape[2]+1 != len(hhl):
      raise IndexError('vertical shape of lwc should match height lvls')

    Nx, Ny, Nz = lwc.shape

    lwc = np.round(lwc, decimals=3)
    reff= np.round(reff,decimals=3)
    reff = np.minimum(np.maximum(reff, reff_min),reff_max)

    dlwc  = xr.DataArray(np.swapaxes(lwc , 0, 1), dims=('ny','nx','nz'))
    dreff = xr.DataArray(np.swapaxes(reff, 0, 1), dims=('ny','nx','nz'))

    D = xr.Dataset({
      'lwc' : xr.DataArray(np.swapaxes(lwc , 0, 1), dims=('ny','nx','nz')),
      'reff': xr.DataArray(np.swapaxes(reff, 0, 1), dims=('ny','nx','nz')),
      'z'   : xr.DataArray(hhl, dims=('nz_lev')), },
      attrs={
        'cldproperties': 3,
        'dx': dx,
        'dy': dy,
        })

    try:
        os.mkdir(os.path.dirname(fname))
    except OSError:
        pass

    print('... writing to file: ... ', fname)
    D.to_netcdf(fname)


from pylab import *
Nx, Ny, Nz = 10, 10, 5

hhl = np.linspace(0,2,Nz+1)
lwc = np.zeros((Nx, Ny, Nz))
reff = np.zeros((Nx, Ny, Nz))

lwc[Nx//2, Ny//2, Nz//2] = 1
reff[lwc>0] = 10

dx = dy = .5

create_librad_netcdf_cldfile('test_cloud.nc', lwc, reff, hhl, dx, dy)

#!/bin/python3
import sys
import os
import shutil
import subprocess
import tempfile
from jinja2 import Template
help_str = """
Generate a LibRadtran input with a template and run it.
Input file is a Temporary file and
Output is saved in a temporary folder
"""

homedir = os.environ['HOME']
libRadtran = os.path.join(homedir, 'libRadtran')

datadict = {
        'data_files_path': os.path.join(libRadtran, 'data'),
        'atmosphere_file': os.path.join(libRadtran, 'data', 'atmmod', 'afglus.dat'),
        'source_solar': os.path.join(libRadtran, 'data', 'atlas_plus_modtran'),
        'mol_abs_param': 'kato2',
        'wavelength_index': '1 32',
        'albedo': 0.15,
        'sza': 60,
        'phi0': 0, # 0 is South

        'rte_solver': 'montecarlo',
        'mc_ipa': True,
        #'mc_tenstream': 'solver_3_10',

        #'rte_solver': 'twostrebe',
        #'ipa_3d': True,

        'mc_sample_grid': '10 10 .5 .5',
        'wc_file_3D': 'test_cloud.nc',
        }

inp_template = Template(
"""
data_files_path {{data_files_path}}
atmosphere_file {{atmosphere_file}}
source solar    {{source_solar}} per_nm

{% if output_process %}
output_process  {{output_process}}
{% else %}
output_process  sum
{% endif %}

mol_abs_param   {{mol_abs_param}}
wavelength_index {{wavelength_index}}

albedo          {{albedo}}
sza             {{sza}}
phi0            {{phi0}}

{% if mc_forward_output %}
mc_forward_output {{mc_forward_output}}
{% else %}
mc_forward_output absorption K_per_day
{% endif %}

{% if zout %}
zout            {{zout}}
{% else %}
zout            all_levels
{% endif %}

rte_solver      {{rte_solver}}
{% if ipa_3d %}
ipa_3d
{% endif %}

{% if mc_tenstream %}
mc_tenstream    {{mc_tenstream}}
{% endif %}

{% if mc_ipa %}
mc_ipa
{% endif %}

{% if mc_photons %}
mc_photons      {{mc_photons}}
{% endif %}

mc_sample_grid  {{mc_sample_grid}}

{% if wc_file_3D %}
wc_file 3D      {{wc_file_3D}}
wc_properties   hu interpolate
{% endif %}

{% if ic_file_3D %}
ic_file 3D      {{ic_file_3D}}
ic_properties   yang interpolate
{% endif %}

{% if outdir %}
mc_basename     {{outdir}}/mc
{% endif %}

#quiet
verbose
""")

tempfile.NamedTemporaryFile
with tempfile.NamedTemporaryFile('w') as fh:
    outdir = tempfile.mkdtemp()
    datadict['outdir'] = outdir

    fh.write(inp_template.render(**datadict))
    fh.flush()

    sp = subprocess.call([os.path.join(libRadtran, 'bin', 'uvspec'), '-f', fh.name])

    shutil.rmtree(outdir)