def create_librad_netcdf_cldfile(fname, lwc, reff, hhl, dx, dy, reff_min=2.51, reff_max=55, round_decimals=None):
    """
    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

    if round_decimals:
        lwc = np.round(lwc, decimals=round_decimals)
        reff= np.round(reff,decimals=round_decimals)

    reff = np.minimum(np.maximum(reff, reff_min),reff_max)

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

    D = xr.Dataset({
      'lwc' : xr.DataArray(lwc , dims=('ny','nx','nz')),
      'reff': xr.DataArray(reff, 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)


def reff_from_lwc_and_N_after_ICON(lwc, N, l_liquid=True, zkap=1.1):
    """
    Compute effective radius from liquid water content and droplet density
    from the implementation in ICON/newcld_optics
    see ECHAM5 documentation (Roeckner et al, MPI report 349)

    lwc :: liquid water content [g m-3],  typically .1
    N   :: droplet density      [1 cm-3], typically 100
    ! zkap => breadth parameter e.g. continental=1.143 (Martin et al.), maritime=1.077(Martin et al.)
    """
    import numpy as np
    rho = 1e3 # water density @ 4degC
    zfact = 1e6 * (3e-9 / (4 * np.pi * rho))**(1./3.) # conversion factor

    if l_liquid:
        return zfact * zkap * (lwc / max(1e-60, N) )**(1./3.)
    else:
        return 83.8 * lwc**0.216
    return None




import xarray as xr
import numpy as np

D=xr.open_dataset('20180115-release-2.3.00-RegA_PR1250m-Default_20131211-11557462-CLv_DOM04_ML_0019_remapcon_ts4_subset600x600pix.nc')

lwc = lwc = D.tot_qc_dia.isel(time=0) #.transpose('lat', 'lon', 'height')
reff = reff_from_lwc_and_N_after_ICON(lwc*1e3,100)

# hydrostatic integration to get heights
# For me, it would make more sense to have the pressure variable be defined on the levels instead of as layer quantities.
# anyway, I'll just put one on top of it that has the same delta_p
pres = D.pres.isel(time=0)[:,::-1,:] # assume that is pressure at levels starting at the bottom
dP = pres.data[:,:-1,:] - pres.data[:,1:,:] # vertical pressue derivative
ndP = np.concatenate((dP, dP[:,-1:,:]), axis=1) # vert. pressure derivative with last value appended

rho = D.rho.isel(time=0)[:,::-1,:] # density starting at the bottom
g = 9.81

dz = ndP / rho / g # vertical extent of layers

height = np.concatenate( (np.zeros_like(dz[:,:1,:]), np.cumsum(dz, axis=1)), axis=1 ) # height levels hydrostatically integrated
height_profile = np.average(height, axis=(0,2)) # horizontally average height profile

# Just to check if that makes sense:
# compute the Liquid water path and plot it:
# lwp = (lwc * dz).sum(dim='height')
# clf(); imshow( lwp.sum(dim='height'), vmin=0, vmax=.1 ); cbar=colorbar(); cbar.set_title('Liq. Water Path [kg/m2]')
# Seems Reasonable with clouds being larger than 100g

# Lat/Lon grid has some distortion but for the sake of LibRadtran computations we neglect thos...
# just assume we a mean cell size
cell_density = np.sqrt( (D.lon.max()-D.lon.min()) * (D.lat.max()-D.lat.min()) / D.lon.size )
m_per_deg = 2. * np.pi * 6371229./ 360.
dx = dy = int( cell_density * m_per_deg)

create_librad_netcdf_cldfile('libradtran_cloud.nc',
        lwc.transpose('lat', 'lon', 'height').data * 1e3,
        reff.transpose('lat', 'lon', 'height').data,
        height_profile * 1e-3,
        dx * 1e-3, dy * 1e-3)