diff --git a/core/core.py b/core/core.py
new file mode 100644
index 0000000000000000000000000000000000000000..a76210123ad9f011129fdd43c729361ad2671a0a
--- /dev/null
+++ b/core/core.py
@@ -0,0 +1,255 @@
+def read_lon_lat(fname: str):
+ """
+ Read geographical position for hdf4 file.
+ Adapted from https://gitlab.phaidra.org/climate/icon_2.3.0_ice-mp_rad_vis/-/blob/master/CloudSat/CloudSat_read.py
+
+ Input:
+ - fname: full path of hdf file (string)
+ Output:
+ - data: xarray dataset with longitude and latitude in deg
+ """
+ import numpy as np
+ from pyhdf import HDF, VS
+ import xarray as xr
+ h = HDF.HDF(fname)
+ vs = h.vstart()
+ # latitude
+ latid = vs.attach(vs.find('Latitude'))
+ latid.setfields('Latitude')
+ nrecs, _, _, _, _ = latid.inquire()
+ latitude = np.array(latid.read(nRec=nrecs)).squeeze()
+ latid.detach()
+ # longitude
+ lonid = vs.attach(vs.find('Longitude'))
+ lonid.setfields('Longitude')
+ nrecs, _, _, _, _ = lonid.inquire()
+ longitude = np.array(lonid.read(nRec=nrecs)).squeeze()
+ lonid.detach()
+ # construct xarray dataset following the dimensions and coordinates of 2B-FLXHR-LIDAR
+ y = np.arange(0,np.size(longitude))+0.5
+ da_lon = xr.DataArray(longitude, name="LON", dims=("y"), coords={"y": y})
+ da_lat = xr.DataArray(latitude, name="LAT", dims=("y"), coords={"y": y})
+ ds = xr.merge([da_lon, da_lat])
+ return ds
+
+
+def read_surfaceheightbin(fname: str):
+ """
+ Read surface height bin. Can be used to mask out values near the surface.
+
+ Input:
+ - fname: full path of hdf file (string)
+ Output:
+ - data: xarray dataset with surface height bin
+ """
+ import numpy as np
+ from pyhdf import HDF, VS
+ import xarray as xr
+ h = HDF.HDF(fname)
+ vs = h.vstart()
+ varid = vs.attach(vs.find('SurfaceHeightBin'))
+ varid.setfields('SurfaceHeightBin')
+ nrecs, _, _, _, _ = varid.inquire()
+ var = np.array(varid.read(nRec=nrecs)).squeeze()
+ varid.detach()
+ # return xarray dataset following the dimensions and coordinates of 2B-FLXHR-LIDAR
+ y = np.arange(0,np.size(var))+0.5
+ da_var = xr.DataArray(var, name="SurfaceHeightBin", dims=("y"), coords={"y": y})
+ return da_var.to_dataset()
+
+
+def read_time(fname: str):
+ """
+ Read time.
+
+ Input:
+ - fname: full path of hdf file (string)
+ Output:
+ - data: xarray dataset with time in datetime64 format
+ """
+ import numpy as np
+ from pyhdf import HDF, VS
+ import xarray as xr
+ from astropy.time import Time
+ h = HDF.HDF(fname)
+ vs = h.vstart()
+ # TAI_start, seconds since 00:00:00 1 Jan 1993
+ varid = vs.attach(vs.find('TAI_start'))
+ varid.setfields('TAI_start')
+ nrecs, _, _, _, _ = varid.inquire()
+ TAI_start = np.array(varid.read(nRec=nrecs)).squeeze()
+ varid.detach()
+ # Profile time
+ varid = vs.attach(vs.find('Profile_time'))
+ varid.setfields('Profile_time')
+ nrecs, _, _, _, _ = varid.inquire()
+ Profile_time = np.array(varid.read(nRec=nrecs)).squeeze()
+ varid.detach()
+ # calculate time using astropy and change to datetim64 format
+ tanom = Time(TAI_start + Profile_time, format="unix_tai", scale="utc")
+ tbase = Time('1993-01-01 00:00:00', scale='utc') # TAI time for CloudSat starts not in 1970 but 1993
+ time = Time(tanom.value+tbase.tai.unix_tai, format="unix_tai", scale="utc")
+ time = time.datetime64
+ # construct xarray dataset following the dimensions and coordinates of 2B-FLXHR-LIDAR
+ y = np.arange(0,np.size(time))+0.5
+ ds = xr.DataArray(time, name="TIME", dims=("y"), coords={"y": y}).to_dataset()
+ return ds
+
+
+def read_height(fname: str):
+ """
+ Read 2d height array for hdf4 file.
+
+ Input:
+ - fname: full path of hdf file (string)
+ Output:
+ - height: xarray dataset with height
+ """
+ import numpy as np
+ from pyhdf.SD import SD, SDC
+ import xarray as xr
+ hdf = SD(fname, SDC.READ)
+ height = hdf.select('Height')
+ height = height[:,:]
+ # construct xarray dataset following the dimensions and coordinates of 2B-FLXHR-LIDAR
+ x = np.arange(0,np.size(height[0]))+0.5
+ y = np.arange(0,np.size(height[:,0]))+0.5
+ da_height = xr.DataArray(height, name="HGT", dims=("y", "xfull"), coords={"y": y, "xfull": x})
+ return da_height.to_dataset()
+
+
+def read_toaswin(fname: str):
+ """
+ Read TOA incoming solar radiation.
+
+ Input:
+ - fname: full path of hdf file (string)
+ Output:
+ - toaswin: xarray dataset with toaswin
+ """
+ from pyhdf import HDF, VS
+ import xarray as xr
+ import numpy as np
+ h = HDF.HDF(fname)
+ vs = h.vstart()
+ varid = vs.attach(vs.find('FD_TOA_IncomingSolar'))
+ varid.setfields('FD_TOA_IncomingSolar')
+ nrecs, _, _, _, _ = varid.inquire()
+ toaswin = np.array(varid.read(nRec=nrecs)).squeeze()
+ varid.detach()
+ # set missing value to nan and rescale
+ missval = toaswin.min()
+ toaswin = 0.1*np.where(toaswin==missval, np.nan, toaswin)
+ # construct xarray dataset following the dimensions and coordinates of 2B-FLXHR-LIDAR
+ y = np.arange(0,np.size(toaswin))+0.5
+ da = xr.DataArray(toaswin, name="TOASWIN", dims=("y"), coords={"y": y})
+ return da.to_dataset()
+
+
+def read_fluxes(fname: str):
+ """
+ Reads radiative fluxes from 2B-FLXHR-LIDAR file.
+
+ Input:
+ - fname: full path of hdf file (string)
+ Output:
+ - xarray dataset with vertically-resolved net longwave and shortwave fluxes.
+ """
+ import rioxarray as rxr
+ import xarray as xr
+ import numpy as np
+
+ data = rxr.open_rasterio(fname)
+
+ # all-sky fluxes
+ FD = data[0]["FD"]
+ missval = FD.min().values
+ FD = 0.1*xr.where(FD==missval, np.nan, FD)
+ FU = data[0]["FU"]
+ missval = FU.min().values
+ FU = 0.1*xr.where(FU==missval, np.nan, FU)
+ FN = (FD - FU).rename("FN")
+
+ # clear-sky fluxes
+ FD_NC = data[0]["FD_NC"]
+ missval = FD_NC.min().values
+ FD_NC = 0.1*xr.where(FD_NC==missval, np.nan, FD_NC)
+ FU_NC = data[0]["FU_NC"]
+ missval = FU_NC.min().values
+ FU_NC = 0.1*xr.where(FU_NC==missval, np.nan, FU_NC)
+ FN_NC = (FD_NC - FU_NC).rename("FN_NC")
+
+ # return all fluxes in one xarray dataset
+ return xr.merge([FN, FN_NC])
+
+
+def read_plevthickness(fname: str):
+ """
+ Calculates pressure level thickness from ECMWF-AUX file.
+
+ Input:
+ - fname: full path of hdf file (string)
+ Output:
+ - xarray dataset with pressure level thickness DP in Pa.
+ """
+ import rioxarray as rxr
+ import xarray as xr
+ import numpy as np
+
+ data = rxr.open_rasterio(fname)
+ P = data["Pressure"].squeeze()
+ missval = P.min().values
+ P = xr.where(P==missval, np.nan, P)
+ DP = 0.5*(P.shift(x=-1) - P.shift(x=+1))
+ DP = DP.rename("DP")
+ return DP
+
+
+def compute_heatingrates(data):
+ """
+ Compute all-sky and clear-sky heating rates from fluxes. Shortwave heating is rescale by the daily-mean
+ insolation.
+
+ Input:
+ - data: xarray dataset with net all-sky (FN) and clear-sky fluxes (FN_NC), TOA incoming shortwave (TOASWIN), pressure level thickness (DP)
+ latitude (LAT) and time (TIME)
+ Output:
+ - output: xarray dataset with all-sky (HR) and clear-sky heating rates (HR_NC) in K/day
+ for rays with zero or nan TOA SW IN, the shortwave heating rate will be nan
+ """
+ import numpy as np
+ import xarray as xr
+ import pandas as pd
+ import sys
+ sys.path.append("/jetfs/home/avoigt/cloudsat-calipso-heating-rates/solar/")
+ import insolation
+ cp = 1004;
+ g = 9.81;
+ HR = -g * data["FN"].diff("x", label="lower") / (cp*data["DP"]) * 86400
+ HR_NC = -g * data["FN_NC"].diff("x", label="lower") / (cp*data["DP"]) * 86400
+
+ # rename
+ HR = HR.rename("HR").rename({"x": "xfull"})
+ HR_NC = HR_NC.rename("HR_NC").rename({"x": "xfull"})
+
+ # split into shortwave and longwave heating
+ HRSW = HR.isel(band=0)
+ HRLW = HR.isel(band=1)
+ HRSW_NC = HR_NC.isel(band=0)
+ HRLW_NC = HR_NC.isel(band=1)
+
+ # rescale shortwave heating rates with daily mean insolation
+ SWIN = data["TOASWIN"] # instantaneous shortwave down at TOA
+ # set SWIN below 50 Wm-2 to nan to remove rays that are close to sunset/sunrise
+ SWIN = xr.where(SWIN<50, np.nan, SWIN)
+ DOY = pd.to_datetime(data.TIME).dayofyear
+ SWIN_daily = insolation.daily_insolation(data.LAT, day=DOY)
+
+ HRSW = HRSW * SWIN_daily/SWIN
+ HRSW_NC = HRSW_NC * SWIN_daily/SWIN
+
+ # combine shortwave and longwave again into one dataarray
+ HR[0] = HRSW; HR[1] = HRLW
+ HR_NC[0] = HRSW_NC; HR_NC[1] = HRLW_NC
+
+ return xr.merge([HR, HR_NC])