diff --git a/.qlty/qlty.toml b/.qlty/qlty.toml
index f911a3979..a08286ac5 100644
--- a/.qlty/qlty.toml
+++ b/.qlty/qlty.toml
@@ -156,5 +156,8 @@ name = "trufflehog"
 [[plugin]]
 name = "yamllint"
 
+[[plugin]]
+name = "black"
+
 #[[plugin]]
 #name = "flake8"
diff --git a/avaframe/com1DFA/DFAfunctionsCython.pxd b/avaframe/com1DFA/DFAfunctionsCython.pxd
index a87b26fa1..7987a1417 100644
--- a/avaframe/com1DFA/DFAfunctionsCython.pxd
+++ b/avaframe/com1DFA/DFAfunctionsCython.pxd
@@ -10,7 +10,7 @@ cpdef (double, double) computeEntMassAndForce(double, double, double, double, do
 
 cpdef double computeDetMass(double, double, double, double)
 
-cpdef double computeResForce(double, double, double, double, double, double, int, int)
+cpdef double computeResForce(double, double, double, double, int, int)
 
 cdef (double, double, double) addArtificialViscosity(double, double, double, double, double, double, double, double,
                                                      int, int, double, double, double, double,
diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx
index d3d0ff1c4..78a389c8d 100644
--- a/avaframe/com1DFA/DFAfunctionsCython.pyx
+++ b/avaframe/com1DFA/DFAfunctionsCython.pyx
@@ -77,7 +77,6 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
   cdef double entDefResistance = cfg.getfloat('entDefResistance')
   cdef double rho = cfg.getfloat('rho')
   cdef double rhoEnt = cfg.getfloat('rhoEnt')
-  cdef double hRes = cfg.getfloat('hRes')
   cdef double gravAcc = cfg.getfloat('gravAcc')
   cdef double xsiVoellmy = cfg.getfloat('xsivoellmy')
   cdef double muVoellmy = cfg.getfloat('muvoellmy')
@@ -371,7 +370,7 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
 
       # adding resistance force due to obstacles
       cResCell = cResRaster[indCellY][indCellX]
-      cResPart = computeResForce(hRes, h, areaPart, rho, cResCell, uMag, explicitFriction, resistanceType)
+      cResPart = computeResForce(areaPart, rho, cResCell, uMag, explicitFriction, resistanceType)
       forceFrict[k] = forceFrict[k] - cResPart
 
       uxArray[k] = ux
@@ -502,16 +501,12 @@ cpdef double computeDetMass(double dt, double detCell,
   return dmDet
 
 
-cpdef double computeResForce(double hRes, double h, double areaPart, double rho, double cResCell,
+cpdef double computeResForce(double areaPart, double rho, double cResCell,
                              double uMag, int explicitFriction, int resistanceType):
   """ compute force component due to resistance
 
   Parameters
   ----------
-  hRes: float
-      resistance height
-  h : float
-      particle flow thickness
   areaPart : float
       particle area
   rho : float
@@ -530,23 +525,15 @@ cpdef double computeResForce(double hRes, double h, double areaPart, double rho,
   cResPart : float
       resistance component for particle
   """
-  cdef double hResEff = hRes
+
   cdef double cRecResPart
-  if(h < hRes):
-      hResEff = h
   # explicit formulation (explicitFriction == 1)
   if explicitFriction == 1:
     if resistanceType == 1:
-      # cRes
-      cRecResPart = - rho * areaPart * hResEff * cResCell * uMag * uMag
-    elif resistanceType == 2:
       # cResH
       cRecResPart = - rho * areaPart * cResCell * uMag * uMag
   elif explicitFriction == 0:
     if resistanceType == 1:
-      # cRes
-      cRecResPart = - rho * areaPart * hResEff * cResCell * uMag
-    elif resistanceType == 2:
       # cResH
       cRecResPart = - rho * areaPart * cResCell * uMag
   return cRecResPart
diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py
index f4cdeef44..4ae908567 100644
--- a/avaframe/com1DFA/com1DFA.py
+++ b/avaframe/com1DFA/com1DFA.py
@@ -1,5 +1,5 @@
 """
-    Main functions for python DFA kernel
+Main functions for python DFA kernel
 """
 
 import copy
@@ -189,7 +189,13 @@ def com1DFAMain(cfgMain, cfgInfo=""):
 
         # postprocessing: writing report, creating plots
         dem, plotDict, reportDictList, simDFNew = com1DFAPostprocess(
-            simDF, tCPUDF, simDFExisting, cfgMain, dem, reportDictList, exportData=exportFlag
+            simDF,
+            tCPUDF,
+            simDFExisting,
+            cfgMain,
+            dem,
+            reportDictList,
+            exportData=exportFlag,
         )
 
         return dem, plotDict, reportDictList, simDFNew
@@ -321,7 +327,12 @@ def com1DFAPostprocess(simDF, tCPUDF, simDFExisting, cfgMain, dem, reportDictLis
         # write report
         reportDictList = gR.checkAndCleanReportDictOnWinIssue872(reportDictList)
 
-        gR.writeReport(reportDir, reportDictList, cfgMain["FLAGS"].getboolean("reportOneFile"), plotDict)
+        gR.writeReport(
+            reportDir,
+            reportDictList,
+            cfgMain["FLAGS"].getboolean("reportOneFile"),
+            plotDict,
+        )
 
     return dem, plotDict, reportDictList, simDFNew
 
@@ -410,9 +421,9 @@ def com1DFACore(cfg, avaDir, cuSimName, inputSimFiles, outDir, simHash=""):
         reportDict["contours"] = contourDictXY
 
     # write text file to Outputs/com1DFA/configurationFilesDone to indicate that this simulation has been performed
-    configFileName = ("%s.ini" % cuSimName)
-    for saveDir in ['configurationFilesDone', 'configurationFilesLatest']:
-        configDir = pathlib.Path(avaDir, "Outputs", "com1DFA", 'configurationFiles', saveDir)
+    configFileName = "%s.ini" % cuSimName
+    for saveDir in ["configurationFilesDone", "configurationFilesLatest"]:
+        configDir = pathlib.Path(avaDir, "Outputs", "com1DFA", "configurationFiles", saveDir)
         with open((configDir / configFileName), "w") as fi:
             fi.write("see directory configurationFiles for info on config")
         fi.close()
@@ -590,7 +601,10 @@ def prepareInputData(inputSimFiles, cfg):
             secondaryReleaseLine["type"] = "Secondary release"
             # check for holes in secondary release area polygons
             gI.checkForMultiplePartsShpArea(
-                cfg["GENERAL"]["avalancheDir"], secondaryReleaseLine, "com1DFA", type="secondary release"
+                cfg["GENERAL"]["avalancheDir"],
+                secondaryReleaseLine,
+                "com1DFA",
+                type="secondary release",
             )
         else:
             message = "No secondary release file found"
@@ -770,7 +784,11 @@ def createReportDict(avaDir, logName, relName, inputSimLines, cfg, reportAreaInf
             "xi": cfgGen["xsivoellmy"],
         }
     elif cfgGen["frictModel"].lower() == "coulomb":
-        reportST["Friction model"] = {"type": "columns", "model": "Coulomb", "mu": cfgGen["mucoulomb"]}
+        reportST["Friction model"] = {
+            "type": "columns",
+            "model": "Coulomb",
+            "mu": cfgGen["mucoulomb"],
+        }
     elif cfgGen["frictModel"].lower() == "wetsnow":
         reportST["Friction model"] = {
             "type": "columns",
@@ -782,9 +800,8 @@ def createReportDict(avaDir, logName, relName, inputSimLines, cfg, reportAreaInf
         reportST["Friction model"] = {
             "type": "columns",
             "model": "spatialVoellmy",
-            "mu file": inputSimLines['muFile'].name,
-            "xi file": inputSimLines['xiFile'].name,
-
+            "mu file": inputSimLines["muFile"].name,
+            "xi file": inputSimLines["xiFile"].name,
         }
 
     # check if secondary release area
@@ -816,7 +833,14 @@ def createReportDict(avaDir, logName, relName, inputSimLines, cfg, reportAreaInf
             }
         )
     if resInfo == "Yes":
-        reportST.update({"Resistance area": {"type": "columns", "Resistance area scenario": resistanceArea}})
+        reportST.update(
+            {
+                "Resistance area": {
+                    "type": "columns",
+                    "Resistance area scenario": resistanceArea,
+                }
+            }
+        )
 
     reportST["Release Area"].update(reportAreaInfo["Release area info"])
 
@@ -994,7 +1018,12 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
         # if no release thickness field or function - set release according to shapefile or ini file
         # this is a list of release rasters that we want to combine
         releaseLine = geoTrans.prepareArea(
-            releaseLine, dem, np.sqrt(2), thList=releaseLine["thickness"], combine=True, checkOverlap=False
+            releaseLine,
+            dem,
+            np.sqrt(2),
+            thList=releaseLine["thickness"],
+            combine=True,
+            checkOverlap=False,
         )
     else:
         # if relTh provided - set release thickness with field or function
@@ -1023,7 +1052,12 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
     releaseLine["header"] = dem["originalHeader"]
     inputSimLines["releaseLine"]["header"] = dem["originalHeader"]
     particles = initializeParticles(
-        cfgGen, releaseLine, dem, inputSimLines=inputSimLines, logName=logName, relThField=relThField
+        cfgGen,
+        releaseLine,
+        dem,
+        inputSimLines=inputSimLines,
+        logName=logName,
+        relThField=relThField,
     )
     particles, fields = initializeFields(cfg, dem, particles, releaseLine)
 
@@ -1040,9 +1074,9 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
     damLine = damCom1DFA.initializeWallLines(cfgGen, dem, damLine, "")
     dem["damLine"] = damLine
     if inputSimLines["damLine"] is None:
-        reportAreaInfo['dam'] = 'No'
+        reportAreaInfo["dam"] = "No"
     else:
-        reportAreaInfo['dam'] = 'Yes'
+        reportAreaInfo["dam"] = "Yes"
 
     # perform initialisation step for redistributing particles
     if cfg["GENERAL"].getboolean("iniStep"):
@@ -1066,7 +1100,12 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
     # get info of simType and whether or not to initialize resistance and entrainment
     simTypeActual = cfgGen["simTypeActual"]
     entrMassRaster, entrEnthRaster, reportAreaInfo = initializeMassEnt(
-        dem, simTypeActual, inputSimLines["entLine"], reportAreaInfo, thresholdPointInPoly, cfgGen
+        dem,
+        simTypeActual,
+        inputSimLines["entLine"],
+        reportAreaInfo,
+        thresholdPointInPoly,
+        cfgGen,
     )
 
     # check if entrainment and release overlap
@@ -1076,10 +1115,18 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
     if secondaryReleaseInfo["flagSecondaryRelease"] == "Yes":
         for secRelRaster in secondaryReleaseInfo["rasterData"]:
             entrMassRaster = geoTrans.checkOverlap(
-                entrMassRaster, secRelRaster, "Entrainment", "Secondary release ", crop=True
+                entrMassRaster,
+                secRelRaster,
+                "Entrainment",
+                "Secondary release ",
+                crop=True,
             )
             entrEnthRaster = geoTrans.checkOverlap(
-                entrEnthRaster, secRelRaster, "Entrainment", "Secondary release ", crop=True
+                entrEnthRaster,
+                secRelRaster,
+                "Entrainment",
+                "Secondary release ",
+                crop=True,
             )
     # surfacic entrainment mass available (unit kg/m²)
     fields["entrMassRaster"] = entrMassRaster
@@ -1089,17 +1136,32 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
 
     log.debug("Initializing resistance area")
     cResRaster, detRaster, reportAreaInfo = initializeResistance(
-        cfgGen, dem, simTypeActual, inputSimLines["resLine"], reportAreaInfo, thresholdPointInPoly
+        cfgGen,
+        dem,
+        simTypeActual,
+        inputSimLines["resLine"],
+        reportAreaInfo,
+        thresholdPointInPoly,
     )
     fields["cResRaster"] = cResRaster
     fields["detRaster"] = detRaster
+    fields["cResRasterOrig"] = cResRaster
+    fields["detRasterOrig"] = detRaster
 
-    for fric in ['mu', 'xi']:
-        if (inputSimLines[fric+'File'] == None) or (cfg['GENERAL']['frictModel'].lower() != 'spatialvoellmy'):
-            fields[fric+'Field'] = np.asarray([[np.nan],[np.nan]])
+    for fric in ["mu", "xi"]:
+        if (inputSimLines[fric + "File"] == None) or (
+            cfg["GENERAL"]["frictModel"].lower() != "spatialvoellmy"
+        ):
+            fields[fric + "Field"] = np.asarray([[np.nan], [np.nan]])
         else:
-            fricField = IOf.readRaster(pathlib.Path(cfg['GENERAL']['avalancheDir'], 'Inputs', cfg['INPUT']['%sFile' % fric]))
-            fields[fric+'Field'] = fricField['rasterData']
+            fricField = IOf.readRaster(
+                pathlib.Path(
+                    cfg["GENERAL"]["avalancheDir"],
+                    "Inputs",
+                    cfg["INPUT"]["%sFile" % fric],
+                )
+            )
+            fields[fric + "Field"] = fricField["rasterData"]
 
     # plot release area scenario
     outCom1DFA.plotReleaseScenarioView(
@@ -1219,7 +1281,16 @@ def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", rel
             hCell = relRaster[indRely, indRelx]
             aCell = areaRaster[indRely, indRelx]
             xPart, yPart, mPart, n, aPart = particleTools.placeParticles(
-                hCell, aCell, indRelx, indRely, csz, massPerPart, nPPK, rng, cfg, ratioArea
+                hCell,
+                aCell,
+                indRelx,
+                indRely,
+                csz,
+                massPerPart,
+                nPPK,
+                rng,
+                cfg,
+                ratioArea,
             )
             nPart = nPart + n
             partPerCell[indRely, indRelx] = n
@@ -1255,6 +1326,7 @@ def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", rel
 
     particles["massPerPart"] = massPerPart
     particles["mTot"] = np.sum(particles["m"])
+    particles["tPlot"] = 0
     particles["h"] = hPartArray
     particles["ux"] = np.zeros(np.shape(hPartArray))
     particles["uy"] = np.zeros(np.shape(hPartArray))
@@ -1292,7 +1364,7 @@ def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", rel
     # where nID is the number of already used IDs
     # (enable tracking of particles even if particles are added or removed)
     # unique identifier for each particle
-    particles["ID"] = np.arange(particles["nPart"],dtype = np.int64)
+    particles["ID"] = np.arange(particles["nPart"], dtype=np.int64)
     # keep track of the identifier (usefull to add identifier to newparticles)
     particles["nID"] = particles["nPart"]
     # keep track of parents (usefull for new particles created after splitting)
@@ -1514,7 +1586,7 @@ def initializeSecRelease(inputSimLines, dem, relRaster, reportAreaInfo):
             np.sqrt(2),
             thList=secondaryReleaseInfo["thickness"],
             combine=False,
-            checkOverlap=False
+            checkOverlap=False,
         )
         # remove overlaping parts of the secondary release area with the main release areas
         noOverlaprasterList = []
@@ -1523,7 +1595,11 @@ def initializeSecRelease(inputSimLines, dem, relRaster, reportAreaInfo):
         ):
             noOverlaprasterList.append(
                 geoTrans.checkOverlap(
-                    secRelRatser, relRaster, "Secondary release " + secRelName, "Release", crop=True
+                    secRelRatser,
+                    relRaster,
+                    "Secondary release " + secRelName,
+                    "Release",
+                    crop=True,
                 )
             )
 
@@ -1587,7 +1663,9 @@ def initializeMassEnt(dem, simTypeActual, entLine, reportAreaInfo, thresholdPoin
         # entrEnthRaster = np.where(tempRaster < 0, tempRaster*cfgGen.getfloat('cpIce'),
         #                           tempRaster*cfgGen.getfloat('cpWtr')/cfgGen.getfloat('hFusion'))
         entrEnthRaster = np.where(
-            entrMassRaster > 0, cfgGen.getfloat("entTempRef") * cfgGen.getfloat("cpIce"), 0
+            entrMassRaster > 0,
+            cfgGen.getfloat("entTempRef") * cfgGen.getfloat("cpIce"),
+            0,
         )
         reportAreaInfo["entrainment"] = "Yes"
     else:
@@ -1628,19 +1706,21 @@ def initializeResistance(cfg, dem, simTypeActual, resLine, reportAreaInfo, thres
     """
     K = cfg.getfloat("detK")
     detrainment = cfg.getboolean("detrainment")
-    detWithoutRes = cfg.getboolean("detWithoutRes")
 
-    ResModel = cfg["ResistanceModel"].lower()
-    if ResModel == "cres":
-        cRes = cfg.getfloat("cRes")
-    if ResModel == "cresh":
-        cRes = cfg.getfloat("cResH")
     # read dem header
     header = dem["originalHeader"]
     ncols = header["ncols"]
     nrows = header["nrows"]
 
     if simTypeActual in ["entres", "res"]:
+        ResModel = cfg["ResistanceModel"].lower()
+        if ResModel == "default":
+            cRes = cfg.getfloat("cResH")
+        else:
+            message = "Resistance model %s not a valid option" % ResModel
+            log.error(message)
+            raise AssertionError(message)
+
         resistanceArea = resLine["fileName"]
         log.info("Initializing resistance area: %s" % (resistanceArea))
         log.info("Resistance area features: %s" % (resLine["Name"]))
@@ -1655,10 +1735,6 @@ def initializeResistance(cfg, dem, simTypeActual, resLine, reportAreaInfo, thres
             log.info("Detrainment (Resistance) area features: %s" % (resLine["Name"]))
             detRaster = K * mask
             reportAreaInfo["detrainment"] = "Yes"
-
-            if detWithoutRes:
-                cResRaster = np.zeros((nrows, ncols))
-                reportAreaInfo["resistance"] = "No"
         else:
             detRaster = np.zeros((nrows, ncols))
             reportAreaInfo["detrainment"] = "No"
@@ -1750,8 +1826,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
     # turn resistance model into integer
     ResModel = cfgGen["ResistanceModel"].lower()
     ResModelsList = [
-        "cres",
-        "cresh",
+        "default",
     ]
     resistanceType = ResModelsList.index(ResModel) + 1
     log.debug("Resistance Model used: %s, %s" % (ResModelsList[resistanceType - 1], resistanceType))
@@ -1762,6 +1837,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
     massEntrained = []
     massDetrained = []
     massTotal = []
+    pfvTimeMax = []
 
     # setup a result fields info data frame to save max values of fields and avalanche front
     resultsDF = setupresultsDF(resTypesLast, cfg["VISUALISATION"].getboolean("createRangeTimeDiagram"))
@@ -1784,7 +1860,10 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
     # export particles properties for visulation
     if cfg["VISUALISATION"].getboolean("writePartToCSV"):
         particleTools.savePartToCsv(
-            cfg["VISUALISATION"]["visuParticleProperties"], [particles], outDir, countParticleCsv=countParticleCsv
+            cfg["VISUALISATION"]["visuParticleProperties"],
+            [particles],
+            outDir,
+            countParticleCsv=countParticleCsv,
         )
         countParticleCsv = countParticleCsv + 1
 
@@ -1842,6 +1921,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
         massDetrained.append(particles["massDetrained"])
         massTotal.append(particles["mTot"])
         timeM.append(t)
+        pfvTimeMax.append(np.nanmax(fields["FV"]))
         # print progress to terminal
         print("time step t = %f s\r" % t, end="")
 
@@ -1886,7 +1966,10 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
             # export particles properties for visulation
             if cfg["VISUALISATION"].getboolean("writePartToCSV"):
                 particleTools.savePartToCsv(
-                    cfg["VISUALISATION"]["visuParticleProperties"], [particles], outDir, countParticleCsv=countParticleCsv
+                    cfg["VISUALISATION"]["visuParticleProperties"],
+                    [particles],
+                    outDir,
+                    countParticleCsv=countParticleCsv,
                 )
                 countParticleCsv = countParticleCsv + 1
 
@@ -1955,6 +2038,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
         "entrained mass": np.sum(massEntrained),
         "detrained mass": np.sum(massDetrained),
         "entrained volume": (np.sum(massEntrained) / cfgGen.getfloat("rhoEnt")),
+        "pfvTimeMax": pfvTimeMax,
     }
 
     # determine if stop criterion is reached or end time
@@ -1987,7 +2071,13 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
         lastTimeStep = t - dt
         # first append final time step
         mtiInfo, dtRangeTime = dtAna.fetchRangeTimeInfo(
-            cfgRangeTime, cfg, dtRangeTime, lastTimeStep, demRT["header"], fields, mtiInfo
+            cfgRangeTime,
+            cfg,
+            dtRangeTime,
+            lastTimeStep,
+            demRT["header"],
+            fields,
+            mtiInfo,
         )
         dtAna.exportData(mtiInfo, cfgRangeTime, "com1DFA")
 
@@ -1997,18 +2087,24 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
 
     if cfg["EXPORTS"].getboolean("exportData"):
         exportFields(cfg, t, fields, dem, outDir, cuSimName, TSave="final")
-        
+
         # export particles dictionaries of saving time steps
         savePartToPickle(particles, outDirData, cuSimName)
     else:
         # fetch contourline info
         contourDictXY = outCom1DFA.fetchContCoors(
-            dem["header"], fields[cfg["VISUALISATION"]["contourResType"]], cfg["VISUALISATION"], cuSimName
+            dem["header"],
+            fields[cfg["VISUALISATION"]["contourResType"]],
+            cfg["VISUALISATION"],
+            cuSimName,
         )
 
     # save contour line for each sim
     contourDictXY = outCom1DFA.fetchContCoors(
-        dem["header"], fields[cfg["VISUALISATION"]["contourResType"]], cfg["VISUALISATION"], cuSimName
+        dem["header"],
+        fields[cfg["VISUALISATION"]["contourResType"]],
+        cfg["VISUALISATION"],
+        cuSimName,
     )
     outDirDataCont = outDir / "contours"
     fU.makeADir(outDirDataCont)
@@ -2017,7 +2113,10 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
     # export particles properties for visulation
     if cfg["VISUALISATION"].getboolean("writePartToCSV"):
         particleTools.savePartToCsv(
-            cfg["VISUALISATION"]["visuParticleProperties"], [particles], outDir, countParticleCsv=countParticleCsv
+            cfg["VISUALISATION"]["visuParticleProperties"],
+            [particles],
+            outDir,
+            countParticleCsv=countParticleCsv,
         )
         countParticleCsv = countParticleCsv + 1
 
@@ -2049,7 +2148,7 @@ def setupresultsDF(resTypes, cfgRangeTime):
 
     resDict = {"timeStep": [0.0]}
     for resT in resTypes:
-        if resT != "particles":
+        if resT != "particles" and resT != "FTDet":
             resDict["max" + resT] = [0.0]
     if cfgRangeTime:
         resDict["rangeList"] = [0.0]
@@ -2083,7 +2182,7 @@ def addMaxValuesToDF(resultsDF, fields, timeStep, resTypes, rangeValue=""):
 
     newLine = []
     for resT in resTypes:
-        if resT != "particles":
+        if resT != "particles" and resT != "FTDet":
             newLine.append(np.nanmax(fields[resT]))
 
     if rangeValue != "":
@@ -2169,6 +2268,12 @@ def writeMBFile(infoDict, avaDir, logName):
     massEntrained = infoDict["massEntrained"]
     massDetrained = infoDict["massDetrained"]
     massTotal = infoDict["massTotal"]
+    massDetrainedTotal = np.zeros(len(massDetrained))
+    for m in range(1, len(massDetrained)):
+        massDetrainedTotal[m] = massDetrainedTotal[m - 1] + massDetrained[m]
+
+    # create mass plot
+    outCom1DFA.massPlot(infoDict, massDetrainedTotal, t, avaDir, logName)
 
     # write mass balance info to log file
     massDir = pathlib.Path(avaDir, "Outputs", "com1DFA")
@@ -2177,8 +2282,14 @@ def writeMBFile(infoDict, avaDir, logName):
         mFile.write("time, current, entrained, detrained\n")
         for m in range(len(t)):
             mFile.write(
-                "%.02f,    %.06f,    %.06f,    %.06f\n"
-                % (t[m], massTotal[m], massEntrained[m], massDetrained[m])
+                "%.02f,    %.06f,    %.06f,   %.06f,    %.06f\n"
+                % (
+                    t[m],
+                    massTotal[m],
+                    massEntrained[m],
+                    massDetrained[m],
+                    massDetrainedTotal[m],
+                )
             )
 
 
@@ -2213,6 +2324,18 @@ def computeEulerTimeStep(cfg, particles, fields, zPartArray0, dem, tCPU, frictTy
     tCPU : dict
         computation time dictionary
     """
+
+    # update cRes and detK rasters according to thresholds of FV and FT
+    particles["tPlot"] = particles["tPlot"] + 1
+    # only if entres or res sim and detrainment is used
+    if cfg["simTypeActual"] in ["entres", "res"] and (
+        cfg["ResistanceModel"].lower() == "default" and cfg.getboolean("detrainment")
+    ):
+        # update resistance area fields using thresholds
+        fields = com1DFATools.updateResCoeffFields(fields, cfg, float(particles["t"]), dem)
+        if debugPlot:
+            outCom1DFA.plotResFields(fields, cfg, particles["tPlot"], dem, particles["mTot"])
+
     # get forces
     startTime = time.time()
     # loop version of the compute force
@@ -2408,7 +2531,10 @@ def trackParticles(cfgTrackPart, dem, particlesList):
     radius = cfgTrackPart.getfloat("radius")
     centerList = cfgTrackPart["centerTrackPartPoint"]
     centerList = centerList.split("|")
-    centerTrackPartPoint = {"x": np.array([float(centerList[0])]), "y": np.array([float(centerList[1])])}
+    centerTrackPartPoint = {
+        "x": np.array([float(centerList[0])]),
+        "y": np.array([float(centerList[1])]),
+    }
     centerTrackPartPoint["x"] = centerTrackPartPoint["x"] - dem["originalHeader"]["xllcenter"]
     centerTrackPartPoint["y"] = centerTrackPartPoint["y"] - dem["originalHeader"]["yllcenter"]
 
@@ -2430,7 +2556,15 @@ def trackParticles(cfgTrackPart, dem, particlesList):
     return particlesList, trackedPartProp, track
 
 
-def readFields(inDir, resType, simName="", flagAvaDir=True, comModule="com1DFA", timeStep="", atol=1.0e-6):
+def readFields(
+    inDir,
+    resType,
+    simName="",
+    flagAvaDir=True,
+    comModule="com1DFA",
+    timeStep="",
+    atol=1.0e-6,
+):
     """Read ascii files within a directory and return List of dictionaries
 
     Parameters
@@ -2505,7 +2639,16 @@ def readFields(inDir, resType, simName="", flagAvaDir=True, comModule="com1DFA",
     return fieldsList, fieldHeader, timeList
 
 
-def exportFields(cfg, timeStep, fields, dem, outDir, cuSimName, TSave="intermediate", resTypesForced=[]):
+def exportFields(
+    cfg,
+    timeStep,
+    fields,
+    dem,
+    outDir,
+    cuSimName,
+    TSave="intermediate",
+    resTypesForced=[],
+):
     """export result fields to Outputs directory according to result parameters and time step
     that can be specified in the configuration file
     option intermediate or final, if final also plotFields for report are exported if intermediate only resTypes are exported
@@ -2541,7 +2684,7 @@ def exportFields(cfg, timeStep, fields, dem, outDir, cuSimName, TSave="intermedi
     if "particles" in resTypesReport:
         resTypesReport.remove("particles")
 
-    if TSave == "final" or TSave == 'initial':
+    if TSave == "final" or TSave == "initial":
         # for last time step we need to add the report fields
         resTypes = list(set(resTypesGen + resTypesReport))
     else:
@@ -2550,7 +2693,11 @@ def exportFields(cfg, timeStep, fields, dem, outDir, cuSimName, TSave="intermedi
     if resTypesForced != []:
         resTypes = resTypesForced
     for resType in resTypes:
-        resField = fields[resType]
+        if resType == "FTDet":
+            dmDet = fields["dmDet"]
+            resField = dmDet / (cfg["GENERAL"].getfloat("rho") * dem["areaRaster"])
+        else:
+            resField = fields[resType]
         if resType == "ppr":
             # convert from Pa to kPa
             resField = resField * 0.001
@@ -2631,7 +2778,10 @@ def prepareVarSimDict(standardCfg, inputSimFiles, variationDict, simNameExisting
         cfgSim = cfgUtils.convertConfigParserToDict(standardCfg)
         # create release scenario name for simulation name
         rel, cfgSim = gI.fetchReleaseFile(
-            inputSimFiles, row._asdict()["releaseScenario"], cfgSim, variationDict["releaseScenario"]
+            inputSimFiles,
+            row._asdict()["releaseScenario"],
+            cfgSim,
+            variationDict["releaseScenario"],
         )
         relName = rel.stem
         if "_" in relName:
@@ -2683,25 +2833,36 @@ def prepareVarSimDict(standardCfg, inputSimFiles, variationDict, simNameExisting
         # check if DEM in Inputs has desired mesh size
         pathToDem = dP.checkRasterMeshSize(cfgSim, inputSimFiles["demFile"], "DEM")
         cfgSim["INPUT"]["DEM"] = pathToDem
-        if cfgSim["GENERAL"]["relThFromFile"] == "True" or cfgSim["GENERAL"]["frictModel"].lower() == "spatialvoellmy":
-            dem = IOf.readRaster(pathlib.Path(cfgSim['GENERAL']['avalancheDir'], 'Inputs', pathToDem))
+        if (
+            cfgSim["GENERAL"]["relThFromFile"] == "True"
+            or cfgSim["GENERAL"]["frictModel"].lower() == "spatialvoellmy"
+        ):
+            dem = IOf.readRaster(pathlib.Path(cfgSim["GENERAL"]["avalancheDir"], "Inputs", pathToDem))
 
         # check if RELTH in Inputs has desired mesh size
         if cfgSim["GENERAL"]["relThFromFile"] == "True":
-            pathToRelTh = dP.checkExtentAndCellSize(cfgSim, inputSimFiles['relThFile'], dem, 'RELTH')
+            pathToRelTh = dP.checkExtentAndCellSize(cfgSim, inputSimFiles["relThFile"], dem, "RELTH")
             cfgSim["INPUT"]["relThFile"] = pathToRelTh
         else:
             cfgSim["INPUT"]["relThFile"] = ""
 
         # check if spatialVoellmy is chosen that friction fields have correct extent
         if cfgSim["GENERAL"]["frictModel"].lower() == "spatialvoellmy":
-            for fric in ['mu', 'xi']:
-                pathToFric = dP.checkExtentAndCellSize(cfgSim, inputSimFiles['%sFile' % fric], dem, fric)
-                cfgSim['INPUT']['%sFile' % fric] = pathToFric
+            for fric in ["mu", "xi"]:
+                pathToFric = dP.checkExtentAndCellSize(cfgSim, inputSimFiles["%sFile" % fric], dem, fric)
+                cfgSim["INPUT"]["%sFile" % fric] = pathToFric
 
         # add info about dam file path to the cfg
-        if cfgSim['GENERAL']['dam'] == 'True' and inputSimFiles['damFile'] != None:
-            cfgSim['INPUT']['DAM'] = str(pathlib.Path('DAM', inputSimFiles['damFile'].name))
+        if cfgSim["GENERAL"]["dam"] == "True" and inputSimFiles["damFile"] is not None:
+            cfgSim["INPUT"]["DAM"] = str(pathlib.Path("DAM", inputSimFiles["damFile"].name))
+
+        # add info about entrainment file path to the cfg
+        if "ent" in row._asdict()["simTypeList"] and inputSimFiles["entFile"] is not None:
+            cfgSim["INPUT"]["entrainment"] = str(pathlib.Path("ENT", inputSimFiles["entFile"].name))
+
+        # add info about entrainment file path to the cfg
+        if "res" in row._asdict()["simTypeList"] and inputSimFiles["resFile"] is not None:
+            cfgSim["INPUT"]["resistance"] = str(pathlib.Path("RES", inputSimFiles["resFile"].name))
 
         # add thickness values if read from shp and not varied
         cfgSim = dP.appendShpThickness(cfgSim)
@@ -2755,7 +2916,10 @@ def prepareVarSimDict(standardCfg, inputSimFiles, variationDict, simNameExisting
             }
             # write configuration file
             cfgUtils.writeCfgFile(
-                cfgSimObject["GENERAL"]["avalancheDir"], com1DFA, cfgSimObject, fileName=simName
+                cfgSimObject["GENERAL"]["avalancheDir"],
+                com1DFA,
+                cfgSimObject,
+                fileName=simName,
             )
         else:
             log.warning("Simulation %s already exists, not repeating it" % simName)
@@ -3004,7 +3168,12 @@ def initializeRelVol(cfg, demVol, releaseFile, radius, releaseType="primary"):
         releaseLine = setThickness(cfg, releaseLine, typeTh)
         # when creating raster from polygon apply release thickness
         releaseLine = geoTrans.prepareArea(
-            releaseLine, demVol, radius, thList=releaseLine["thickness"], combine=True, checkOverlap=False
+            releaseLine,
+            demVol,
+            radius,
+            thList=releaseLine["thickness"],
+            combine=True,
+            checkOverlap=False,
         )
 
         # compute release volume using raster and dem area
@@ -3014,16 +3183,16 @@ def initializeRelVol(cfg, demVol, releaseFile, radius, releaseType="primary"):
 
 
 def saveContToPickle(contourDictXY, outDir, cuSimName):
-    """ save contourline x, y coordinates dictionary to a pickle
-
-        Parameters
-        ------------
-        contourDictXY: dict
-            dictionary with key simName and dict with x, y coordinates of contour line of specified level
-        outDir: pathlib path
-            path to dir where pickle shall be saved
-        cuSimName: str
-            name of current simulation where this contourline is derived from
+    """save contourline x, y coordinates dictionary to a pickle
+
+    Parameters
+    ------------
+    contourDictXY: dict
+        dictionary with key simName and dict with x, y coordinates of contour line of specified level
+    outDir: pathlib path
+        path to dir where pickle shall be saved
+    cuSimName: str
+        name of current simulation where this contourline is derived from
     """
 
     fi = open(outDir / ("contourDictXY_%s.pickle" % (cuSimName)), "wb")
diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini
index 91e191039..257a8b06e 100644
--- a/avaframe/com1DFA/com1DFACfg.ini
+++ b/avaframe/com1DFA/com1DFACfg.ini
@@ -353,25 +353,31 @@ entTempRef = -10.0
 # J/kgK (ice)
 cpIce = 2050.0
 
-#++++++++++++ Resistance force parameters
-# Choose one resistance model:
-# cRes : F = - cRes *hEff * rho * A* u²; hEff = min(hRes,h) [cRes] = 1/m
-# cResH : F = - cResH * rho * A* u², [cResH] =[]
-ResistanceModel = cRes
-# height of the obstacles [m]
-hRes = 10
-# resistance parameters for different models
-cRes = 0.003
-cResH = 0.003
-
-#++++++++++++ Detrainment parameter
-# Flag if snow detrainment at obstacles (in resistance area) is computed
-# detrainment can be added to resistance force (DetAndRes = True) or computed solely (DetWithoutRes = True)
-detrainment = False
-detWithoutRes = False
+#++++++++++++ Resistance model
+# default setup:
+ResistanceModel = default
+# At each time step, ResistanceModel default applies increased friction and optional detrainment (see below). Only relevant in resistance areas.
+# NOTE: development setup; parameter values need more testing and calibration!!
+
+#  parameter for increased friction in resistance areas
+cResH = 0.01
+
+# Apply detrainment in resistance areas in accordance with the flow thickness and flow velocity thresholds specified below.
+# if False - only increased friction is applied in resistance areas
+detrainment = True
+
 # detrainment parameter defined by Feistl et al. (2014)
-# value need to be refined
-detK = 20  
+detK = 5
+
+# thresholds if detrainment is set to True
+# FV OR FT below min thresholds: apply only detrainment. no increased friction
+# FV AND FT within min and max thresholds: no detrainment, only apply increased friction
+# FV OR FT above max thresholds: no detrainment and no increased friction applied
+forestVMin = 6.
+forestThMin = 0.6
+forestVMax = 40.
+forestThMax = 10.
+
 
 #++++++++++++ Entrainment Erosion Energy
 # Used to determine speed loss via energy loss due to entrained mass
diff --git a/avaframe/com1DFA/com1DFATools.py b/avaframe/com1DFA/com1DFATools.py
index 068996166..a837fbbf6 100644
--- a/avaframe/com1DFA/com1DFATools.py
+++ b/avaframe/com1DFA/com1DFATools.py
@@ -1,5 +1,5 @@
 """
-    Tools specific to the com1DFA computational kernel
+Tools specific to the com1DFA computational kernel
 """
 
 import configparser
@@ -8,6 +8,7 @@
 import logging
 import math
 import pathlib
+import numpy as np
 
 from deepdiff import DeepDiff
 
@@ -17,7 +18,8 @@
 from avaframe.com1DFA import com1DFA
 from avaframe.in1Data import getInput as gI
 from avaframe.in3Utils import cfgUtils
-from avaframe.in3Utils import fileHandlerUtils as fU
+from avaframe.in2Trans import rasterUtils as IOf
+
 
 # create local logger
 # change log level in calling module to DEBUG to see log messages
@@ -328,7 +330,8 @@ def initializeInputs(avalancheDir, cleanRemeshedRasters):
     # if need to reproduce exactly the hash - need to be strings with exactly the same number of digits!!
     # searchCfgFiles=True enables to search for preformed sims if run has been interrupted
     simDFExisting, simNameExisting = cfgUtils.readConfigurationInfoFromDone(
-        avalancheDir, specDir='',
+        avalancheDir,
+        specDir="",
     )
 
     # fetch input data - dem, release-, entrainment- and resistance areas (and secondary release areas)
@@ -380,3 +383,69 @@ def checkCfgInfoType(cfgInfo):
         raise AssertionError(message)
 
     return typeCfgInfo
+
+
+def updateResCoeffFields(fields, cfg, t, dem):
+    """update fields of cRes and detK, coefficients of resistance parameter and detrainment parameter
+    according to the thresholds of FV and FT
+    if FV OR FT below min thresholds -> apply detrainment in that area
+    if FV AND FT above min thresholds AND below max thresholds -> apply cResH (additional friction) in that area
+    if FV OR FT above max thresholds -> no impact of resistance area
+
+    Parameters
+    ------------
+    fields: dict
+        dictionary with cResRasterOrig, cResRaster and detRasterOrig, detRaster fields
+    cfg: configparser object
+        configuration of com1DFA, thresholds
+    t: float
+        time step
+    dem: dict
+        dictionary with info on DEM
+
+    Returns
+    ------------
+    fields: dict
+        updated cRes and detK fields, required to write raster with header
+    """
+
+    # fetch cRes and detK raster and thresholds for FV and FT
+    cResOrig = fields["cResRasterOrig"].copy()
+    detOrig = fields["detRasterOrig"].copy()
+    vMin = cfg.getfloat("forestVMin")
+    thMin = cfg.getfloat("forestThMin")
+    vMax = cfg.getfloat("forestVMax")
+    thMax = cfg.getfloat("forestThMax")
+
+    # create new rasters using FV, FT and thresholds to mask
+    detRasterInt = np.where(((fields["FV"] <= vMin) | (fields["FT"] <= thMin)), detOrig, 0.0)
+    detRaster = np.where(((fields["FV"] > vMax) | (fields["FT"] > thMax)), 0, detRasterInt)
+    cResRasterInt = np.where(((fields["FV"] > vMin) & (fields["FT"] > thMin)), cResOrig, 0.0)
+    cResRaster = np.where(((fields["FV"] > vMax) | (fields["FT"] > thMax)), 0, cResRasterInt)
+
+    if len(np.where((detRaster > 0) & (cResRaster > 0))[0]) > 0:
+        message = "Detrainment and increased friction within same cell!"
+        log.error(message)
+        raise AssertionError(message)
+
+    # if max thresholds are exceeded: forest destroyed remove forest
+    lTh = len(np.where((fields["FV"] > vMax) | (fields["FT"] > thMax))[0])
+    if lTh > 0:
+        cResOrig = np.where(((fields["FV"] > vMax) | (fields["FT"] > thMax)), 0, cResOrig)
+        detOrig = np.where(((fields["FV"] > vMax) | (fields["FT"] > thMax)), 0, detOrig)
+        fields["cResRasterOrig"] = cResOrig
+        fields["detRasterOrig"] = detOrig
+        log.info(
+            "Resistance area removed %d cells because FV or FT exceeded %.2f ms-1, %.2f m"
+            % (lTh, vMax, thMax)
+        )
+        outFileRes = pathlib.Path(
+            cfg["avalancheDir"], "Outputs", "com1DFA", "reports", ("resArea_t%.2f" % t)
+        )
+        IOf.writeResultToRaster(dem["originalHeader"], fields["cResRasterOrig"], outFileRes, flip=True)
+
+    # update fields dictionary
+    fields["cResRaster"] = cResRaster
+    fields["detRaster"] = detRaster
+
+    return fields
diff --git a/avaframe/com1DFA/deriveParameterSet.py b/avaframe/com1DFA/deriveParameterSet.py
index a25e43e02..2fbaefc53 100644
--- a/avaframe/com1DFA/deriveParameterSet.py
+++ b/avaframe/com1DFA/deriveParameterSet.py
@@ -1,5 +1,5 @@
 """
-    Create dictionary for parameter variations
+Create dictionary for parameter variations
 """
 
 import logging
@@ -149,6 +149,7 @@ def checkResType(fullCfg, section, key, value):
             "TA",
             "particles",
             "dmDet",
+            "FTDet",
         ]
         message = "The parameter % s is not a valid resType. It will not be saved"
         newResType = []
@@ -385,15 +386,27 @@ def checkThicknessSettings(cfg, thName):
     thRV = thName + "RangeVariation"
     thPV = thName + "PercentVariation"
     thRCiV = thName + "RangeFromCiVariation"
-    flagsList = [cfg["GENERAL"][thRV] != "", cfg["GENERAL"][thPV] != "", cfg["GENERAL"][thRCiV] != ""]
+    flagsList = [
+        cfg["GENERAL"][thRV] != "",
+        cfg["GENERAL"][thPV] != "",
+        cfg["GENERAL"][thRCiV] != "",
+    ]
 
     if sum(flagsList) > 1:
-        message = "Only one variation type is allowed - check %s and %s, %s" % (thRV, thPV, thRCiV)
+        message = "Only one variation type is allowed - check %s and %s, %s" % (
+            thRV,
+            thPV,
+            thRCiV,
+        )
         log.error(message)
         raise AssertionError(message)
 
     if cfg["GENERAL"].getboolean(thFile) and (sum(flagsList) > 0):
-        message = "RelThFromFile is True - no variation allowed: check %s, %s or %s" % (thRV, thPV, thRCiV)
+        message = "RelThFromFile is True - no variation allowed: check %s, %s or %s" % (
+            thRV,
+            thPV,
+            thRCiV,
+        )
         log.error(message)
         raise AssertionError(message)
 
@@ -907,7 +920,15 @@ def checkExtentAndCellSize(cfg, inputFile, dem, fileType):
             "Field %s interpolated onto DEM extent and corresponding spatial resolution, "
             "cellSize changed from %.2f to %.2f; difference of llcenter was in x: %.2f, in y: %.2f m"
             "and urcenter was in x: %.2f, in y %.2f"
-            % (fileType, cellSizeOld, dem["header"]["cellsize"], diffX0, diffY0, diffX1, diffY1)
+            % (
+                fileType,
+                cellSizeOld,
+                dem["header"]["cellsize"],
+                diffX0,
+                diffY0,
+                diffX1,
+                diffY1,
+            )
         )
 
         # save to inputField
diff --git a/avaframe/com5SnowSlide/com5SnowSlideCfg.ini b/avaframe/com5SnowSlide/com5SnowSlideCfg.ini
index a455c7c8a..24482d8b2 100644
--- a/avaframe/com5SnowSlide/com5SnowSlideCfg.ini
+++ b/avaframe/com5SnowSlide/com5SnowSlideCfg.ini
@@ -61,7 +61,5 @@ minDistCohesion = 1.0e-3
 cohesiveSurfaceTension = 50000
 
 #++++++++++++ Resistance force parameters
-# height of obstacles [m]
-hRes = 100
 # characteristic coefficient
-cRes = 30000
+cResH = 30000
diff --git a/avaframe/out1Peak/outPlotAllPeak.py b/avaframe/out1Peak/outPlotAllPeak.py
index 65dcfe4e5..724a3e87b 100644
--- a/avaframe/out1Peak/outPlotAllPeak.py
+++ b/avaframe/out1Peak/outPlotAllPeak.py
@@ -50,7 +50,7 @@ def plotAllPeakFields(avaDir, cfgFLAGS, modName, demData=""):
     # Load all infos on simulations
     avaDir = pathlib.Path(avaDir)
     inputDir = avaDir / "Outputs" / modName / "peakFiles"
-    inDir = avaDir / 'Inputs'
+    inDir = avaDir / "Inputs"
     peakFilesDF = fU.makeSimDF(inputDir, avaDir=avaDir)
 
     if demData == "":
@@ -60,7 +60,9 @@ def plotAllPeakFields(avaDir, cfgFLAGS, modName, demData=""):
     else:
         # check if nodata_value is found and if replace with nans for plotting
         demDataField = np.where(
-            demData["rasterData"] == demData["header"]["nodata_value"], np.nan, demData["rasterData"]
+            demData["rasterData"] == demData["header"]["nodata_value"],
+            np.nan,
+            demData["rasterData"],
         )
     demField = demDataField
 
@@ -79,12 +81,11 @@ def plotAllPeakFields(avaDir, cfgFLAGS, modName, demData=""):
 
     # Loop through peakFiles and generate plot
     for m in range(len(peakFilesDF["names"])):
-
         # Load names and paths of peakFiles
         name = peakFilesDF["names"][m]
         fileName = peakFilesDF["files"][m]
         resType = peakFilesDF["resType"][m]
-        simType = peakFilesDF['simType'][m]
+        simType = peakFilesDF["simType"][m]
 
         log.debug("now plot %s:" % (fileName))
 
@@ -94,7 +95,6 @@ def plotAllPeakFields(avaDir, cfgFLAGS, modName, demData=""):
         # make sure to remove the Outputs folder if you want to regenerate the plot
         # this enables to append simulations to an already existing output without regenerating all plots
         if not plotName.is_file():
-
             # Figure  shows the result parameter data
             fig, ax = plt.subplots(figsize=(pU.figW, pU.figH))
 
@@ -102,43 +102,68 @@ def plotAllPeakFields(avaDir, cfgFLAGS, modName, demData=""):
             cellSize = peakFilesDF["cellSize"][m]
 
             # add peak field data now
-            ax, rowsMinPlot, colsMinPlot, extentCellCorners = addConstrainedDataField(fileName, resType, demField, ax, cellSize)
+            ax, rowsMinPlot, colsMinPlot, extentCellCorners = addConstrainedDataField(
+                fileName, resType, demField, ax, cellSize
+            )
 
             # Add background map
             if cfgFLAGS.getboolean("showOnlineBackground"):
                 rasterInfo = IOf.readRaster(fileName, noDataToNan=True)
                 providers = ctx.providers.flatten()
-                srcCrs = rasterInfo['header']['crs']
+                srcCrs = rasterInfo["header"]["crs"]
                 if srcCrs is None:
-                    message = 'chosen basemap: %s not applicable for CRS: %s' % (str(cfgFLAGS["mapProvider"]), srcCrs)
+                    message = "chosen basemap: %s not applicable for CRS: %s" % (
+                        str(cfgFLAGS["mapProvider"]),
+                        srcCrs,
+                    )
                     log.error(message)
                     raise AssertionError(message)
-                ctx.add_basemap(ax, crs=srcCrs, source=providers[str(cfgFLAGS["mapProvider"])], zorder=2)
+                ctx.add_basemap(
+                    ax,
+                    crs=srcCrs,
+                    source=providers[str(cfgFLAGS["mapProvider"])],
+                    zorder=2,
+                )
 
             # if entrainment or resistance area is considered in simulation, show extent of entrainment or resistance area
-            colorOutline = {'ent': 'white', 'res': 'green'}
-            for sType in ['ent', 'res']:
+            colorOutline = {"ent": "white", "res": "green"}
+            for sType in ["ent", "res"]:
                 if sType in simType:
                     sFile, sInfo = gI.getAndCheckInputFiles(inDir, sType.upper(), sType, fileExt="shp")
-                    if sInfo != 'No':
+                    if sInfo != "No":
                         sarea = gpd.read_file(sFile)
-                        sarea.plot(ax=ax, zorder=12, edgecolor=colorOutline[sType], linewidth=2, facecolor="none",
-                                       label=('%s area' % sType), alpha=0.8)
+                        sarea.plot(
+                            ax=ax,
+                            zorder=12,
+                            edgecolor=colorOutline[sType],
+                            linewidth=2,
+                            facecolor="none",
+                            label=("%s area" % sType),
+                            alpha=0.8,
+                        )
 
             # set limit to axis from constrainedData
             ax.set_xlim(extentCellCorners[0], extentCellCorners[1])
             ax.set_ylim(extentCellCorners[2], extentCellCorners[3])
 
             # if available zoom into area provided by crop shp file in Inputs/CROPSHAPE
-            cropFile, cropInfo = gI.getAndCheckInputFiles(inDir, "POLYGONS", "cropFile", fileExt="cropshape.shp")
-            if cropInfo != 'No':
+            cropFile, cropInfo = gI.getAndCheckInputFiles(
+                inDir, "POLYGONS", "cropFile", fileExt="shp", fileSuffix="_cropshape"
+            )
+            if cropInfo != "No":
                 focus = gpd.read_file(cropFile)
-                focus.plot(ax=ax, zorder=12, edgecolor="red", linewidth=2, facecolor="none", alpha=0)
+                focus.plot(
+                    ax=ax,
+                    zorder=12,
+                    edgecolor="red",
+                    linewidth=2,
+                    facecolor="none",
+                    alpha=0,
+                )
                 extent = focus.total_bounds
                 ax.set_xlim(extent[0], extent[2])
                 ax.set_ylim(extent[1], extent[3])
 
-
             # add title, labels and ava Info
             title = str("%s" % name)
             ax.set_title(title)
@@ -153,32 +178,32 @@ def plotAllPeakFields(avaDir, cfgFLAGS, modName, demData=""):
     return plotDict
 
 
-def addConstrainedDataField(fileName, resType, demField, ax, cellSize, alpha=1., oneColor=''):
-    """ find fileName data, constrain data and demField to where there is data,
-        create colormap, define extent, add hillshade contours, add to axes
-        and add colorbar
-
-        Parameters
-        -----------
-        fileName: pathlib path
-            path to data
-        resType: str
-            name of result variable type
-        demField: numpy ndarray
-            array of dem data
-        ax: matplotlib axes object
-            axes where to add data plot to
-        cellSize: float
-            cellSize of data
-        alpha: float
-            from 0 transparent to 1 opaque for plot of constrained data
-        oneColor: str
-            optional to add a color for a single color for field
-
-        Return
-        --------
-        ax: matplotlib axes object
-            axes updated
+def addConstrainedDataField(fileName, resType, demField, ax, cellSize, alpha=1.0, oneColor=""):
+    """find fileName data, constrain data and demField to where there is data,
+    create colormap, define extent, add hillshade contours, add to axes
+    and add colorbar
+
+    Parameters
+    -----------
+    fileName: pathlib path
+        path to data
+    resType: str
+        name of result variable type
+    demField: numpy ndarray
+        array of dem data
+    ax: matplotlib axes object
+        axes where to add data plot to
+    cellSize: float
+        cellSize of data
+    alpha: float
+        from 0 transparent to 1 opaque for plot of constrained data
+    oneColor: str
+        optional to add a color for a single color for field
+
+    Return
+    --------
+    ax: matplotlib axes object
+        axes updated
     """
 
     # Load data
@@ -187,8 +212,8 @@ def addConstrainedDataField(fileName, resType, demField, ax, cellSize, alpha=1.,
 
     # constrain data to where there is data
     rowsMin, rowsMax, colsMin, colsMax = pU.constrainPlotsToData(data, cellSize)
-    dataConstrained = data[rowsMin:rowsMax + 1, colsMin:colsMax + 1]
-    demConstrained = demField[rowsMin:rowsMax + 1, colsMin:colsMax + 1]
+    dataConstrained = data[rowsMin : rowsMax + 1, colsMin : colsMax + 1]
+    demConstrained = demField[rowsMin : rowsMax + 1, colsMin : colsMax + 1]
 
     data = np.ma.masked_where(dataConstrained == 0.0, dataConstrained)
     dataConstrained = np.ma.masked_where(dataConstrained == 0.0, dataConstrained)
@@ -200,32 +225,50 @@ def addConstrainedDataField(fileName, resType, demField, ax, cellSize, alpha=1.,
 
     # choose colormap
     cmap, col, ticks, norm = pU.makeColorMap(
-        pU.colorMaps[resType], np.amin(data), np.amax(data), continuous=pU.contCmap)
+        pU.colorMaps[resType], np.amin(data), np.amax(data), continuous=pU.contCmap
+    )
     cmap.set_bad(alpha=0)
     # uncomment this to set the under value for discrete cmap transparent
     # cmap.set_under(alpha=0)
 
     # set extent in meters using cellSize and llcenter location
-    extentCellCenters, extentCellCorners, rowsMinPlot, rowsMaxPlot, colsMinPlot, colsMaxPlot= pU.createExtent(rowsMin, rowsMax, colsMin, colsMax, raster['header'])
+    (
+        extentCellCenters,
+        extentCellCorners,
+        rowsMinPlot,
+        rowsMaxPlot,
+        colsMinPlot,
+        colsMaxPlot,
+    ) = pU.createExtent(rowsMin, rowsMax, colsMin, colsMax, raster["header"])
 
     # add DEM hillshade with contour lines
     _, _ = pU.addHillShadeContours(ax, demConstrained, cellSize, extentCellCenters)
 
     # add peak field data
-    if oneColor != '':
-        dataOneColor = np.where(dataConstrained > 0.0, np.amax(data)*0.25, np.nan)
-        ax.imshow(dataOneColor, cmap=oneColor, norm=norm, extent=extentCellCorners, origin="lower", aspect="equal", zorder=4,
-                        alpha=alpha)
+    if oneColor != "":
+        dataOneColor = np.where(dataConstrained > 0.0, np.amax(data) * 0.25, np.nan)
+        ax.imshow(
+            dataOneColor,
+            cmap=oneColor,
+            norm=norm,
+            extent=extentCellCorners,
+            origin="lower",
+            aspect="equal",
+            zorder=4,
+            alpha=alpha,
+        )
     else:
-        im1 = ax.imshow(dataConstrained,
-                        cmap=cmap,
-                        norm=norm,
-                        extent=extentCellCorners,
-                        origin="lower",
-                        aspect="equal",
-                        zorder=4,
-                        alpha=alpha)
-        if len(np.nonzero(data)[0]) > 0.:
+        im1 = ax.imshow(
+            dataConstrained,
+            cmap=cmap,
+            norm=norm,
+            extent=extentCellCorners,
+            origin="lower",
+            aspect="equal",
+            zorder=4,
+            alpha=alpha,
+        )
+        if len(np.nonzero(data)[0]) > 0.0:
             # add Colorbar
             fig = ax.get_figure()
             divider = make_axes_locatable(ax)
@@ -263,7 +306,6 @@ def plotAllFields(avaDir, inputDir, outDir, unit="", constrainData=True):
 
     # Loop through peakFiles and generate plot
     for filename in peakFiles:
-
         # Load data
         raster = IOf.readRaster(filename)
         data = raster["rasterData"]
@@ -291,10 +333,12 @@ def plotAllFields(avaDir, inputDir, outDir, unit="", constrainData=True):
 
         if constrainData:
             rowsMin, rowsMax, colsMin, colsMax = pU.constrainPlotsToData(data, cellSize)
-            dataConstrained = data[rowsMin:rowsMax + 1, colsMin:colsMax + 1]
+            dataConstrained = data[rowsMin : rowsMax + 1, colsMin : colsMax + 1]
             data = np.ma.masked_where(dataConstrained == 0.0, dataConstrained)
             # set extent in meters using cellSize and llcenter location
-            extentCellCenters, extentCellCorners, _, _, _, _ = pU.createExtent(rowsMin, rowsMax, colsMin, colsMax, header)
+            extentCellCenters, extentCellCorners, _, _, _, _ = pU.createExtent(
+                rowsMin, rowsMax, colsMin, colsMax, header
+            )
             im1 = ax.imshow(
                 data,
                 cmap=cmap,
@@ -306,7 +350,12 @@ def plotAllFields(avaDir, inputDir, outDir, unit="", constrainData=True):
         else:
             extentCellCenters, extentCellCorners = pU.createExtentMinMax(data, header, originLLCenter=True)
             im1 = ax.imshow(
-                data, cmap=cmap, norm=norm, extent=extentCellCorners, origin="lower", aspect=nx / ny
+                data,
+                cmap=cmap,
+                norm=norm,
+                extent=extentCellCorners,
+                origin="lower",
+                aspect=nx / ny,
             )
 
         pU.addColorBar(im1, ax, ticks, unit)
diff --git a/avaframe/out3Plot/outCom1DFA.py b/avaframe/out3Plot/outCom1DFA.py
index ef97d4c3f..8e690cdf0 100644
--- a/avaframe/out3Plot/outCom1DFA.py
+++ b/avaframe/out3Plot/outCom1DFA.py
@@ -16,13 +16,13 @@
 import avaframe.in3Utils.fileHandlerUtils as fU
 
 cfgMain = cfgUtils.getGeneralConfig()
-cfgFlags = cfgMain['FLAGS']
+cfgFlags = cfgMain["FLAGS"]
 # create local logger
 log = logging.getLogger(__name__)
 
 
 def plotTrackParticle(outDirData, particlesList, trackedPartProp, cfg, dem, cuSimName):
-    """ Plot time series of tracked particles
+    """Plot time series of tracked particles
 
     Parameters
     ----------
@@ -41,59 +41,60 @@ def plotTrackParticle(outDirData, particlesList, trackedPartProp, cfg, dem, cuSi
         name of current simulation
 
     """
-    cfgTrackPart = cfg['TRACKPARTICLES']
-    radius = cfgTrackPart.getfloat('radius')
-    centerList = cfgTrackPart['centerTrackPartPoint']
-    centerList = centerList.split('|')
-    center = {'x': np.array([float(centerList[0])]),
-              'y': np.array([float(centerList[1])])}
-    center, _ = geoTrans.projectOnRaster(dem, center, interp='bilinear')
-    time = trackedPartProp['t']
+    cfgTrackPart = cfg["TRACKPARTICLES"]
+    radius = cfgTrackPart.getfloat("radius")
+    centerList = cfgTrackPart["centerTrackPartPoint"]
+    centerList = centerList.split("|")
+    center = {
+        "x": np.array([float(centerList[0])]),
+        "y": np.array([float(centerList[1])]),
+    }
+    center, _ = geoTrans.projectOnRaster(dem, center, interp="bilinear")
+    time = trackedPartProp["t"]
 
     # do some ploting
-    fig = plt.figure(figsize=(pU.figW*3, pU.figH*2))
-    fig.suptitle('Tracked particles')
+    fig = plt.figure(figsize=(pU.figW * 3, pU.figH * 2))
+    fig.suptitle("Tracked particles")
     ax1 = plt.subplot(221)
-    ax1 = addDem2Plot(ax1, dem, what='slope')
-    circle1 = plt.Circle((center['x'], center['y']), radius, color='r')
-    ax1.plot(trackedPartProp['x'], trackedPartProp['y'])
+    ax1 = addDem2Plot(ax1, dem, what="slope")
+    circle1 = plt.Circle((center["x"], center["y"]), radius, color="r")
+    ax1.plot(trackedPartProp["x"], trackedPartProp["y"])
     ax1.add_patch(circle1)
-    ax1.set_xlabel('x [m]')
-    ax1.set_ylabel('y [m]')
-    ax1.set_title('Trajectory')
+    ax1.set_xlabel("x [m]")
+    ax1.set_ylabel("y [m]")
+    ax1.set_title("Trajectory")
 
     ax2 = plt.subplot(222)
-    ax2.plot(time, trackedPartProp['m'])
-    ax2.set_xlabel('t [s]')
-    ax2.set_ylabel('m [kg]')
-    ax2.set_title('Mass')
+    ax2.plot(time, trackedPartProp["m"])
+    ax2.set_xlabel("t [s]")
+    ax2.set_ylabel("m [kg]")
+    ax2.set_title("Mass")
 
     ax3 = plt.subplot(223)
-    velocity = DFAtls.norm(trackedPartProp['ux'], trackedPartProp['uy'],
-                           trackedPartProp['uz'])
+    velocity = DFAtls.norm(trackedPartProp["ux"], trackedPartProp["uy"], trackedPartProp["uz"])
     ax3.plot(time, velocity)
-    ax3.set_xlabel('t [s]')
-    ax3.set_ylabel('v [m/s]')
-    ax3.set_title('Velocity')
+    ax3.set_xlabel("t [s]")
+    ax3.set_ylabel("v [m/s]")
+    ax3.set_title("Velocity")
 
     ax4 = plt.subplot(224)
-    ax4.plot(time, trackedPartProp['h'])
-    ax4.set_xlabel('t [s]')
-    ax4.set_ylabel('h [m]')
-    ax4.set_title('Flow thickness')
+    ax4.plot(time, trackedPartProp["h"])
+    ax4.set_xlabel("t [s]")
+    ax4.set_ylabel("h [m]")
+    ax4.set_title("Flow thickness")
 
     pathDict = {}
-    pathDict['pathResult'] = outDirData
-    outFileName = 'trackedParticles_%s' % cuSimName
+    pathDict["pathResult"] = outDirData
+    outFileName = "trackedParticles_%s" % cuSimName
     pU.saveAndOrPlot(pathDict, outFileName, fig)
 
-    if cfgFlags.getboolean('showPlot'):
+    if cfgFlags.getboolean("showPlot"):
         fig2 = plt.figure()
         ax1 = plt.subplot(111)
         for count in range(len(particlesList)):
             particles = particlesList[count]
             ax1 = updateTrackPart(particles, ax1, dem)
-            ax1 = addDem2Plot(ax1, dem, what='slope')
+            ax1 = addDem2Plot(ax1, dem, what="slope")
             plt.pause(0.1)
         plt.show()
 
@@ -107,7 +108,7 @@ def plotTrackParticle(outDirData, particlesList, trackedPartProp, cfg, dem, cuSi
 
 
 def plotTrackParticleAcceleration(outDirData, trackedPartProp, cfg, cuSimName):
-    """ Plot time series of tracked particles
+    """Plot time series of tracked particles
     Parameters
     ----------
     outDirData: str
@@ -124,77 +125,76 @@ def plotTrackParticleAcceleration(outDirData, trackedPartProp, cfg, cuSimName):
 
     # do some ploting
     fig = plt.figure(figsize=(pU.figW, pU.figH))
-    fig.suptitle('Tracked particles acceleration')
+    fig.suptitle("Tracked particles acceleration")
     ax1 = plt.subplot(111)
 
-    ax1.plot(trackedPartProp['t'], trackedPartProp['uAcc'])
-    ax1.set_xlabel('time step [s]')
-    ax1.set_ylabel('acceleration [ms-2]')
+    ax1.plot(trackedPartProp["t"], trackedPartProp["uAcc"])
+    ax1.set_xlabel("time step [s]")
+    ax1.set_ylabel("acceleration [ms-2]")
 
     pathDict = {}
-    pathDict['pathResult'] = outDirData
-    outFileName = 'trackedParticles_acceleration_ %s' % cuSimName
+    pathDict["pathResult"] = outDirData
+    outFileName = "trackedParticles_acceleration_ %s" % cuSimName
     pU.saveAndOrPlot(pathDict, outFileName, fig)
 
 
 def plotAllPartAcc(outDirData, particlesList, cfg, Tsave, cuSimName):
-    """ Plot time series of tracked particles
-        Parameters
-        ----------
-        outDirData: str
-            path to output directory
-        particlesList: list
-            list with dictionary of all particles time steps
-            particles
-        cfg : dict
-            configuration read from ini file
-        Tsave: list
-            list of saving time step info
-        cuSimName: str
-            name of current sim
+    """Plot time series of tracked particles
+    Parameters
+    ----------
+    outDirData: str
+        path to output directory
+    particlesList: list
+        list with dictionary of all particles time steps
+        particles
+    cfg : dict
+        configuration read from ini file
+    Tsave: list
+        list of saving time step info
+    cuSimName: str
+        name of current sim
     """
 
     # initialize fig
     fig = plt.figure(figsize=(pU.figW, pU.figH))
-    fig.suptitle('Tracked particles acceleration')
+    fig.suptitle("Tracked particles acceleration")
     ax1 = plt.subplot(111)
 
-    uAcc = np.zeros((len(particlesList), particlesList[0]['nPart']))
-    timeStep = np.asarray([p['t'] for p in particlesList])
-    for idx in particlesList[0]['ID']:
-        uAcc[:, idx] = np.asarray([p['uAcc'][idx] for p in particlesList])
+    uAcc = np.zeros((len(particlesList), particlesList[0]["nPart"]))
+    timeStep = np.asarray([p["t"] for p in particlesList])
+    for idx in particlesList[0]["ID"]:
+        uAcc[:, idx] = np.asarray([p["uAcc"][idx] for p in particlesList])
         ax1.plot(Tsave, uAcc[:, idx])
-    ax1.set_xlabel('time step [s]')
-    ax1.set_ylabel('acceleration [ms-2]')
+    ax1.set_xlabel("time step [s]")
+    ax1.set_ylabel("acceleration [ms-2]")
 
     pathDict = {}
-    pathDict['pathResult'] = outDirData
-    outFileName = 'allparticles_acceleration_%s' % cuSimName
+    pathDict["pathResult"] = outDirData
+    outFileName = "allparticles_acceleration_%s" % cuSimName
     pU.saveAndOrPlot(pathDict, outFileName, fig)
 
 
 def updateTrackPart(particles, ax, dem):
-    """Update axes with particles (tracked particles are highlighted in red)
-    """
+    """Update axes with particles (tracked particles are highlighted in red)"""
 
-    header = dem['header']
-    xllc = header['xllcenter']
-    yllc = header['yllcenter']
+    header = dem["header"]
+    xllc = header["xllcenter"]
+    yllc = header["yllcenter"]
 
-    X = particles['x'] + xllc
-    Y = particles['y'] + yllc
+    X = particles["x"] + xllc
+    Y = particles["y"] + yllc
 
     ax.clear()
-    ax.set_title('t=%.2f s' % particles['t'])
-    variable = particles['trackedParticles']
+    ax.set_title("t=%.2f s" % particles["t"])
+    variable = particles["trackedParticles"]
     # set range and steps of colormap
     cc = np.where(variable == 1, True, False)
-    ax.scatter(X, Y, c='b', cmap=None, marker='.')
-    ax.scatter(X[cc], Y[cc], c='r', cmap=None, marker='.', s=5)
+    ax.scatter(X, Y, c="b", cmap=None, marker=".")
+    ax.scatter(X[cc], Y[cc], c="r", cmap=None, marker=".", s=5)
     return ax
 
 
-def addParticles2Plot(particles, ax, dem, whatS='m', whatC='h', colBarResType=''):
+def addParticles2Plot(particles, ax, dem, whatS="m", whatC="h", colBarResType=""):
     """Update axes with particles
 
     Parameters
@@ -211,33 +211,37 @@ def addParticles2Plot(particles, ax, dem, whatS='m', whatC='h', colBarResType=''
 
     """
 
-    header = dem['header']
-    xllc = header['xllcenter']
-    yllc = header['yllcenter']
+    header = dem["header"]
+    xllc = header["xllcenter"]
+    yllc = header["yllcenter"]
 
-    X = particles['x'] + xllc
-    Y = particles['y'] + yllc
+    X = particles["x"] + xllc
+    Y = particles["y"] + yllc
     cmap = pU.cmapT
-    if colBarResType != '':
-        unit = pU.cfgPlotUtils['unit%s' % colBarResType]
+    if colBarResType != "":
+        unit = pU.cfgPlotUtils["unit%s" % colBarResType]
     else:
-        unit = ''
+        unit = ""
     variableC = particles[whatC]
     variableS = particles[whatS]
-    minMax = np.nanmax(variableS)-np.nanmin(variableS)
+    minMax = np.nanmax(variableS) - np.nanmin(variableS)
     if minMax > 0:
-        variableS = ((variableS-np.nanmin(variableS))/(np.nanmax(variableS)-np.nanmin(variableS)) + 1)*pU.ms
+        variableS = (
+            (variableS - np.nanmin(variableS)) / (np.nanmax(variableS) - np.nanmin(variableS)) + 1
+        ) * pU.ms
     else:
         variableS = pU.ms
-    cmap, _, ticks, norm = pU.makeColorMap(cmap, np.amin(variableC), np.amax(variableC), continuous=pU.contCmap)
+    cmap, _, ticks, norm = pU.makeColorMap(
+        cmap, np.amin(variableC), np.amax(variableC), continuous=pU.contCmap
+    )
     # set range and steps of colormap
-    sc = ax.scatter(X, Y, c=variableC, s=variableS, cmap=cmap, marker='.', zorder=15)
+    sc = ax.scatter(X, Y, c=variableC, s=variableS, cmap=cmap, marker=".", zorder=15)
     cb = pU.addColorBar(sc, ax, ticks, unit)
     return ax, cb
 
 
-def addDem2Plot(ax, dem, what='slope', extent='', origHeader=False):
-    """ Add dem to the background of a plot
+def addDem2Plot(ax, dem, what="slope", extent="", origHeader=False):
+    """Add dem to the background of a plot
 
     Parameters
     ----------
@@ -254,52 +258,73 @@ def addDem2Plot(ax, dem, what='slope', extent='', origHeader=False):
         if True use originalHeader and not header
     """
     if origHeader:
-        header = dem['originalHeader']
+        header = dem["originalHeader"]
     else:
-        header = dem['header']
-    ncols = header['ncols']
-    nrows = header['nrows']
-    xllc = header['xllcenter']
-    yllc = header['yllcenter']
-    csz = header['cellsize']
-    xArray = np.linspace(xllc, xllc+(ncols-1)*csz, ncols)
-    yArray = np.linspace(yllc, yllc+(nrows-1)*csz, nrows)
+        header = dem["header"]
+    ncols = header["ncols"]
+    nrows = header["nrows"]
+    xllc = header["xllcenter"]
+    yllc = header["yllcenter"]
+    csz = header["cellsize"]
+    xArray = np.linspace(xllc, xllc + (ncols - 1) * csz, ncols)
+    yArray = np.linspace(yllc, yllc + (nrows - 1) * csz, nrows)
     cmap = pU.cmapGreys
-    cmap.set_bad(color='white')
+    cmap.set_bad(color="white")
 
-    if what == 'slope':
-        value = dem['Nz'] / DFAtls.norm(dem['Nx'], dem['Ny'], dem['Nz'])
-    elif what == 'z':
-        value = dem['rasterData']
-    elif what == 'hillshade':
+    if what == "slope":
+        value = dem["Nz"] / DFAtls.norm(dem["Nx"], dem["Ny"], dem["Nz"])
+    elif what == "z":
+        value = dem["rasterData"]
+    elif what == "hillshade":
         ls = pU.LightSource(azdeg=pU.azimuthDegree, altdeg=pU.elevationDegree)
-        value = dem['rasterData']
-        value = ls.hillshade(value, vert_exag=pU.vertExag, dx=value.shape[1],
-                             dy=value.shape[0])
+        value = dem["rasterData"]
+        value = ls.hillshade(value, vert_exag=pU.vertExag, dx=value.shape[1], dy=value.shape[0])
     else:
-        value = dem['rasterData']
-    if extent == '':
+        value = dem["rasterData"]
+    if extent == "":
         extent = [xArray.min(), xArray.max(), yArray.min(), yArray.max()]
 
     # here extent is in data coordinates (here cell centers)
-    ref0, im = pU.NonUnifIm(ax, xArray, yArray, value, 'x [m]', 'y [m]',
-                            # extent=[2400, 2700, YY.min(), YY.max()],
-                            extent=extent,
-                            cmap=cmap, norm=None, zorder=0)
-    CS = ax.contour(xArray, yArray, dem['rasterData'], levels=pU.hillshadeContLevs, colors='k', linewidths=0.5,
-                    zorder=3)
+    ref0, im = pU.NonUnifIm(
+        ax,
+        xArray,
+        yArray,
+        value,
+        "x [m]",
+        "y [m]",
+        # extent=[2400, 2700, YY.min(), YY.max()],
+        extent=extent,
+        cmap=cmap,
+        norm=None,
+        zorder=0,
+    )
+    CS = ax.contour(
+        xArray,
+        yArray,
+        dem["rasterData"],
+        levels=pU.hillshadeContLevs,
+        colors="k",
+        linewidths=0.5,
+        zorder=3,
+    )
 
     # add labels
     ax.clabel(CS, CS.levels, inline=True, fontsize=8, zorder=3)
 
     # add info box with indication of label meaning
-    pU.putInfoBox(ax, '- elevation [m]', location='upperLeft', color='gray',
-                  hAlignment='left', alphaF=1.0)
+    pU.putInfoBox(
+        ax,
+        "- elevation [m]",
+        location="upperLeft",
+        color="gray",
+        hAlignment="left",
+        alphaF=1.0,
+    )
     return ax
 
 
 def plotParticles(particlesList, cfg, dem):
-    """ Plot particles on dem
+    """Plot particles on dem
 
     Parameters
     ----------
@@ -311,20 +336,20 @@ def plotParticles(particlesList, cfg, dem):
         dem dictionary with normal information
 
     """
-    if cfgFlags.getboolean('showPlot'):
+    if cfgFlags.getboolean("showPlot"):
         for count in range(len(particlesList)):
             fig2 = plt.figure()
             ax1 = plt.subplot(111)
             ax1.clear()
             particles = particlesList[count]
-            ax1.set_title('t=%.2f s' % particles['t'])
-            ax1 = addParticles2Plot(particles, ax1, dem, whatS='h', whatC='m')
-            ax1 = addDem2Plot(ax1, dem, what='slope')
+            ax1.set_title("t=%.2f s" % particles["t"])
+            ax1 = addParticles2Plot(particles, ax1, dem, whatS="h", whatC="m")
+            ax1 = addDem2Plot(ax1, dem, what="slope")
             plt.show()
 
 
 def addResult2Plot(ax, header, rasterData, resType, colorbar=True, contour=False):
-    """ Add raster data to a plot
+    """Add raster data to a plot
 
     Parameters
     ----------
@@ -338,33 +363,47 @@ def addResult2Plot(ax, header, rasterData, resType, colorbar=True, contour=False
     colorbar: bool
         If true add the colorbar
     """
-    xllc = header['xllcenter']
-    yllc = header['yllcenter']
-    csz = header['cellsize']
-    rowsMin, rowsMax, colsMin, colsMax, rasterData = pU.constrainPlotsToData(rasterData, csz, extentOption=False,
-                                                                             constrainedData=True)
+    xllc = header["xllcenter"]
+    yllc = header["yllcenter"]
+    csz = header["cellsize"]
+    rowsMin, rowsMax, colsMin, colsMax, rasterData = pU.constrainPlotsToData(
+        rasterData, csz, extentOption=False, constrainedData=True
+    )
     # here it is xStart:xEnd where xEnd is colsMax not colsMax-1 as colsMax should still be included where
     # it would be ncols-1 if ncols is the total number of cols
-    xArray = np.linspace(xllc+colsMin*csz, xllc+colsMax*csz, colsMax-colsMin+1)
-    yArray = np.linspace(yllc+rowsMin*csz, yllc+rowsMax*csz, rowsMax-rowsMin+1)
+    xArray = np.linspace(xllc + colsMin * csz, xllc + colsMax * csz, colsMax - colsMin + 1)
+    yArray = np.linspace(yllc + rowsMin * csz, yllc + rowsMax * csz, rowsMax - rowsMin + 1)
     extent = [xArray.min(), xArray.max(), yArray.min(), yArray.max()]
-    unit = pU.cfgPlotUtils['unit%s' % resType]
-    contourLevels = pU.cfgPlotUtils['contourLevels%s' % resType]
+    unit = pU.cfgPlotUtils["unit%s" % resType]
+    contourLevels = pU.cfgPlotUtils["contourLevels%s" % resType]
     contourLevels = fU.splitIniValueToArraySteps(contourLevels)
-    cmap, _, ticks, norm = pU.makeColorMap(pU.colorMaps[resType], np.amin(rasterData), np.amax(rasterData),
-                                           continuous=pU.contCmap)
+    cmap, _, ticks, norm = pU.makeColorMap(
+        pU.colorMaps[resType],
+        np.amin(rasterData),
+        np.amax(rasterData),
+        continuous=pU.contCmap,
+    )
     cmap.set_bad(alpha=0)
     rasterData = np.ma.masked_where(rasterData == 0, rasterData)
 
-    ref0, im = pU.NonUnifIm(ax, xArray, yArray, rasterData, 'x [m]', 'y [m]',
-                            extent=extent,
-                            cmap=cmap, norm=norm, zorder=9)
+    ref0, im = pU.NonUnifIm(
+        ax,
+        xArray,
+        yArray,
+        rasterData,
+        "x [m]",
+        "y [m]",
+        extent=extent,
+        cmap=cmap,
+        norm=norm,
+        zorder=9,
+    )
     if colorbar:
         cb = pU.addColorBar(im, ax, ticks, unit)
     else:
         cb = None
     if contour:
-        CS = ax.contour(xArray, yArray, rasterData, levels=contourLevels, colors='k')
+        CS = ax.contour(xArray, yArray, rasterData, levels=contourLevels, colors="k")
         ax.clabel(CS, inline=1, fontsize=8, zorder=10)
     else:
         CS = None
@@ -372,182 +411,388 @@ def addResult2Plot(ax, header, rasterData, resType, colorbar=True, contour=False
 
 
 def createContourPlot(reportDictList, avalancheDir, simDF):
-    """ create a contour line plot of all simulations of current run
-
-        Parameters
-        -----------
-        reportDictList: list
-            list of com1DFA dictionary with info on contour dicts for each sim
-        avalancheDir: str or pathlib path
-            path to avalanche directory
-        simDF: pandas DataFrame
-            dataframe with one row per simulation performed and its parameter settings
-
-        Returns
-        -------
-        reportDictList: list
-            updated reportDictList - deleted contours dict
+    """create a contour line plot of all simulations of current run
+
+    Parameters
+    -----------
+    reportDictList: list
+        list of com1DFA dictionary with info on contour dicts for each sim
+    avalancheDir: str or pathlib path
+        path to avalanche directory
+    simDF: pandas DataFrame
+        dataframe with one row per simulation performed and its parameter settings
+
+    Returns
+    -------
+    reportDictList: list
+        updated reportDictList - deleted contours dict
     """
 
-    modName = 'com1DFA'
+    modName = "com1DFA"
     contourD = {}
 
     # fetch coordinates of contour line for each sim in reportDict and create contourD
     for cont in reportDictList:
-        contourD[list(cont['contours'].keys())[0]] = cont['contours'][list(cont['contours'].keys())[0]]
-        del cont['contours']
+        contourD[list(cont["contours"].keys())[0]] = cont["contours"][list(cont["contours"].keys())[0]]
+        del cont["contours"]
 
     # create contours directory in Ouputs to save plot
-    contourPlotDir = pathlib.Path(avalancheDir, 'Outputs', modName, 'contours')
+    contourPlotDir = pathlib.Path(avalancheDir, "Outputs", modName, "contours")
     fU.makeADir(contourPlotDir)
 
     # create pathDict and fetch simName of first sim found
-    pathDict = {'pathResult': contourPlotDir, 'plotScenario': pathlib.Path(avalancheDir).stem,
-        'avaDir': avalancheDir}
-    simName1 = simDF['simName'].iloc[0]
+    pathDict = {
+        "pathResult": contourPlotDir,
+        "plotScenario": pathlib.Path(avalancheDir).stem,
+        "avaDir": avalancheDir,
+    }
+    simName1 = simDF["simName"].iloc[0]
     # check if consistent settings throughout all sims
-    if len(np.unique(simDF['contourResType'])) != 1 or len(np.unique(simDF['thresholdValue'])) != 1:
-        log.warning('Contour result type or thresholdValue are not identical for all sims performed - so cannot create contour plot')
+    if len(np.unique(simDF["contourResType"])) != 1 or len(np.unique(simDF["thresholdValue"])) != 1:
+        log.warning(
+            "Contour result type or thresholdValue are not identical for all sims performed - so cannot create contour plot"
+        )
         return reportDictList, False
 
     # generate plot
-    oQ.plotContours(contourD,
-                    contourD[simName1]['contourResType'],
-                    contourD[simName1]['thresholdValue'],
-                    pathDict)
+    oQ.plotContours(
+        contourD,
+        contourD[simName1]["contourResType"],
+        contourD[simName1]["thresholdValue"],
+        pathDict,
+    )
 
-    log.info('Saved contour plot of %d sims to %s' % (len(simDF), contourPlotDir))
-    log.info('Sim names are: %s' % simDF['simName'].to_list())
+    log.info("Saved contour plot of %d sims to %s" % (len(simDF), contourPlotDir))
+    log.info("Sim names are: %s" % simDF["simName"].to_list())
 
     return reportDictList, True
 
 
 def fetchContCoors(demHeader, flowF, cfgVisu, simName):
-    """ fetch coordinates of contour line
-
-        Parameters
-        ------------
-        demHeader: dict
-            dictionary of dem nrows, ncols, cellsize
-        flowF: np.ndarray
-            field data used to compute contour line
-        cfgVisu: configparser object
-            configuration settings for visualisation, here used:
-            contourResType, thresholdValue
-        simName: str
-            simName
-
-        Rertuns
-        --------
-        contDictXY: dict
-            dictionary with simName and subDict with x, y coordinates of contour line
-            contourResType and thresholdValue
+    """fetch coordinates of contour line
+
+    Parameters
+    ------------
+    demHeader: dict
+        dictionary of dem nrows, ncols, cellsize
+    flowF: np.ndarray
+        field data used to compute contour line
+    cfgVisu: configparser object
+        configuration settings for visualisation, here used:
+        contourResType, thresholdValue
+    simName: str
+        simName
+
+    Rertuns
+    --------
+    contDictXY: dict
+        dictionary with simName and subDict with x, y coordinates of contour line
+        contourResType and thresholdValue
     """
     # create coordinate grid first
     xGrid, yGrid, _, _ = geoTrans.makeCoordGridFromHeader(demHeader)
 
     # fetch contour line
-    contourDictXYLines = pU.fetchContourCoords(xGrid, yGrid, flowF, cfgVisu.getfloat('thresholdValue'))
+    contourDictXYLines = pU.fetchContourCoords(xGrid, yGrid, flowF, cfgVisu.getfloat("thresholdValue"))
 
     # setup dict
     contDictXY = {simName: contourDictXYLines}
-    contDictXY[simName]['contourResType'] = cfgVisu['contourResType']
-    contDictXY[simName]['thresholdValue'] = cfgVisu.getfloat('thresholdValue')
+    contDictXY[simName]["contourResType"] = cfgVisu["contourResType"]
+    contDictXY[simName]["thresholdValue"] = cfgVisu.getfloat("thresholdValue")
 
     return contDictXY
 
 
-def plotReleaseScenarioView(avaDir, releaseLine, relThField, reportAreaInfo, dem, titleFig, cuSimName, inputSimLines):
-    """ plot release polygon, area with thickness on dem hillshade
-        saved to avaDir/Outputs/com1DFA/reports
-
-        Parameters
-        ------------
-        avaDir: str
-            path to ava directory
-        dem: dict
-            dict with dem header and data
-        releaseLine, damLine, entLine, resLine, secondaryReleaseLine: dict
-            dict with raster of release, dam, entrainment, resistance, secondary release line and x,y coors
-        reportAreaInfo: dict
-            info with Yes or No for entrainment, resistance, secRelArea, dam
-        titleFig: str
-            title of figure
-        cuSimName: str
-            name of simulation
+def plotReleaseScenarioView(
+    avaDir,
+    releaseLine,
+    relThField,
+    reportAreaInfo,
+    dem,
+    titleFig,
+    cuSimName,
+    inputSimLines,
+):
+    """plot release polygon, area with thickness on dem hillshade
+    saved to avaDir/Outputs/com1DFA/reports
+
+    Parameters
+    ------------
+    avaDir: str
+        path to ava directory
+    dem: dict
+        dict with dem header and data
+    releaseLine, damLine, entLine, resLine, secondaryReleaseLine: dict
+        dict with raster of release, dam, entrainment, resistance, secondary release line and x,y coors
+    reportAreaInfo: dict
+        info with Yes or No for entrainment, resistance, secRelArea, dam
+    titleFig: str
+        title of figure
+    cuSimName: str
+        name of simulation
 
     """
 
-    ny = releaseLine['rasterData'].shape[0]
-    nx = releaseLine['rasterData'].shape[1]
-    Ly = ny * dem['originalHeader']['cellsize']
-    Lx = nx * dem['originalHeader']['cellsize']
-    xL = dem['originalHeader']['xllcenter']
-    yL = dem['originalHeader']['yllcenter']
-    originCells = dem['header']['cellsize'] * 0.5
+    ny = releaseLine["rasterData"].shape[0]
+    nx = releaseLine["rasterData"].shape[1]
+    Ly = ny * dem["originalHeader"]["cellsize"]
+    Lx = nx * dem["originalHeader"]["cellsize"]
+    xL = dem["originalHeader"]["xllcenter"]
+    yL = dem["originalHeader"]["yllcenter"]
+    originCells = dem["header"]["cellsize"] * 0.5
     if len(relThField) == 0:
-        releaseF = releaseLine['rasterData'].copy()
+        releaseF = releaseLine["rasterData"].copy()
     else:
-        releaseF = np.where(releaseLine['rasterData'] > 0, relThField, 0)
+        releaseF = np.where(releaseLine["rasterData"] > 0, relThField, 0)
 
     rField = np.ma.masked_where(releaseF == 0.0, releaseF)
 
     # choose colormap
     cmap1, col, ticks, norm = pU.makeColorMap(
-        pU.colorMaps['pft'], np.amin(releaseLine['rasterData']), np.amax(releaseLine['rasterData']),
-        continuous=pU.contCmap)
-    cmap1.set_bad(alpha=0.)
+        pU.colorMaps["pft"],
+        np.amin(releaseLine["rasterData"]),
+        np.amax(releaseLine["rasterData"]),
+        continuous=pU.contCmap,
+    )
+    cmap1.set_bad(alpha=0.0)
     # extent taking into account the 0.5*cellSize for imshow plot if meters is used
-    extentCells = [xL - originCells, xL + Lx - originCells, yL + Ly - originCells, yL - originCells]
+    extentCells = [
+        xL - originCells,
+        xL + Lx - originCells,
+        yL + Ly - originCells,
+        yL - originCells,
+    ]
     # extent given in cell center coordinates
     extentDem = [xL, xL + Lx, yL + Ly, yL]
 
     # create figure
     fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(pU.figW, pU.figH))
-    addDem2Plot(ax, dem, what='hillshade', extent=extentDem, origHeader=True)
+    addDem2Plot(ax, dem, what="hillshade", extent=extentDem, origHeader=True)
     im1 = ax.imshow(rField, extent=extentCells, cmap=cmap1)
     handles = []
-    relArea = gpd.read_file(inputSimLines['releaseLine']['file'])
+    relArea = gpd.read_file(inputSimLines["releaseLine"]["file"])
     relArea.plot(ax=ax, edgecolor="darkblue", linewidth=2, facecolor="none")
-    relPatch = Patch(color='darkblue', label='release')
+    relPatch = Patch(color="darkblue", label="release")
     handles.append(relPatch)
 
     count = 1
-    if reportAreaInfo['resistance'] == 'Yes':
-        resArea = gpd.read_file(inputSimLines['resLine']['fileName'])
+    if reportAreaInfo["resistance"] == "Yes":
+        resArea = gpd.read_file(inputSimLines["resLine"]["fileName"])
         resArea.plot(ax=ax, edgecolor="green", linewidth=2, facecolor="none")
-        resPatch = Patch(color='green', label='resistance')
+        resPatch = Patch(color="green", label="resistance")
         handles.append(resPatch)
         count = count + 1
-    if reportAreaInfo['entrainment'] == 'Yes':
-        entArea = gpd.read_file(inputSimLines['entLine']['fileName'])
+    if reportAreaInfo["entrainment"] == "Yes":
+        entArea = gpd.read_file(inputSimLines["entLine"]["fileName"])
         entArea.plot(ax=ax, edgecolor="lightblue", linewidth=2, facecolor="none")
-        entPatch = Patch(color='lightblue', label='entrainment')
+        entPatch = Patch(color="lightblue", label="entrainment")
         handles.append(entPatch)
         count = count + 1
-    if reportAreaInfo['secRelArea'] != 'No':
-        secRelArea = gpd.read_file(inputSimLines['secondaryReleaseLine']['fileName'])
+    if reportAreaInfo["secRelArea"] != "No":
+        secRelArea = gpd.read_file(inputSimLines["secondaryReleaseLine"]["fileName"])
         secRelArea.plot(ax=ax, edgecolor="blue", linewidth=2, facecolor="none")
-        secRelPatch = Patch(color='blue', label='secondary release')
+        secRelPatch = Patch(color="blue", label="secondary release")
         handles.append(secRelPatch)
         count = count + 1
-    if reportAreaInfo['dam'] == 'Yes':
-        damArea = gpd.read_file(inputSimLines['damLine']['fileName'][0])
+    if reportAreaInfo["dam"] == "Yes":
+        damArea = gpd.read_file(inputSimLines["damLine"]["fileName"][0])
         damArea.plot(ax=ax, edgecolor="orange", linewidth=2, facecolor="none")
-        damPatch = Patch(color='orange', label='dam')
+        damPatch = Patch(color="orange", label="dam")
         handles.append(damPatch)
         count = count + 1
 
-    ax.set_aspect('equal')
-    cax = ax.inset_axes([1.04, 0.0, 0.05, 1.])
-    pU.addColorBar(im1, ax, ticks, 'm', cax=cax)
-    plt.legend(handles=handles, fontsize=8, loc='upper center', bbox_to_anchor=(0.5, -0.15),
-                 ncol=int(np.ceil(count/2)))
+    ax.set_aspect("equal")
+    cax = ax.inset_axes([1.04, 0.0, 0.05, 1.0])
+    pU.addColorBar(im1, ax, ticks, "m", cax=cax)
+    plt.legend(
+        handles=handles,
+        fontsize=8,
+        loc="upper center",
+        bbox_to_anchor=(0.5, -0.15),
+        ncol=int(np.ceil(count / 2)),
+    )
     plt.title(titleFig, fontsize=8)
     pU.putAvaNameOnPlot(ax, avaDir)
 
     # save and or plot
-    plotName = 'releaseScenario_%s' % cuSimName
-    outDir = pathlib.Path(avaDir, 'Outputs', 'com1DFA', 'reports')
+    plotName = "releaseScenario_%s" % cuSimName
+    outDir = pathlib.Path(avaDir, "Outputs", "com1DFA", "reports")
+    fU.makeADir(outDir)
+    pU.saveAndOrPlot({"pathResult": outDir}, plotName, fig)
+
+
+def plotResFields(fields, cfg, tPlot, dem, totalMass):
+    """Plot active resistance area fields for com1DFArun
+
+    Parameters
+    -----------
+    fields: dict
+        dictionary with fields of FV, FT, cResRaster, detKRaster of current time step
+    cfg: configparser object
+        configuration settings of com1DFA - FV, FT thresholds for resistance processes
+    tPlot: int
+        integer counting for plot naming, ascending order with time
+    dem: dict
+        dictionary with info on DEM
+    """
+
+    # load thresholds
+    vMin = cfg.getfloat("forestVMin")
+    thMin = cfg.getfloat("forestThMin")
+    vMax = cfg.getfloat("forestVMax")
+    thMax = cfg.getfloat("forestThMax")
+
+    # fetch DEM info
+    header = dem["originalHeader"]
+    # create extent for cell centers lower left to upper right in meters
+    extentCellCenters, extentCellCorners = pU.createExtentMinMax(
+        dem["rasterData"], header, originLLCenter=True
+    )
+
+    # create figure
+    fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(pU.figW * 2, pU.figH))
+
+    # mask 0 values in arrays for plots
+    cRes = np.flipud(np.where(fields["cResRaster"] == 0, np.nan, fields["cResRaster"]))
+    det = np.flipud(np.where(fields["detRaster"] == 0, np.nan, fields["detRaster"]))
+    fv = np.flipud(np.where(fields["FV"] == 0, np.nan, fields["FV"]))
+    ft = np.flipud(np.where(fields["FT"] == 0, np.nan, fields["FT"]))
+    dmDet = np.flipud(np.where(fields["dmDet"] == 0, np.nan, fields["dmDet"]))
+    FTDet = dmDet / (cfg.getfloat("rho") * dem["areaRaster"])
+
+    # plot resistance model parameter and detrainment parameter
+    im1 = ax[0, 0].imshow(cRes, vmin=0, vmax=cfg.getfloat("cRes"), extent=extentCellCenters, zorder=4)
+    im2 = ax[0, 1].imshow(det, vmin=0, vmax=cfg.getfloat("detK"), extent=extentCellCenters, zorder=4)
+
+    # set title
+    ax[0, 0].set_title("cResRaster")
+    ax[0, 1].set_title("detRaster")
+    # add colorbar
+    fig.colorbar(im1, ax=ax[0, 0])
+    fig.colorbar(im2, ax=ax[0, 1])
+
+    # fetch colormaps for FV and FT colorcode values below and above thresholds, mask 0 values
+    cmap1, col, ticks, norm = pU.makeColorMap(pU.colorMaps["FV"], vMin, vMax, continuous=pU.contCmap)
+    cmap2, col, ticks, norm = pU.makeColorMap(pU.colorMaps["FV"], thMin, thMax, continuous=pU.contCmap)
+    cmap1.set_under(color="lightblue")
+    cmap1.set_over(color="red")
+    cmap2.set_under(color="lightblue")
+    cmap2.set_over(color="red")
+    cmap1.set_bad(color="k", alpha=0.0)
+    cmap2.set_bad(color="k", alpha=0.0)
+
+    # add FV and FT fields
+    im3 = ax[1, 0].imshow(fv, cmap=cmap1, vmin=vMin, vmax=vMax, extent=extentCellCenters, zorder=4)
+    im4 = ax[1, 1].imshow(ft, cmap=cmap2, vmin=thMin, vmax=thMax, extent=extentCellCenters, zorder=4)
+    # set title
+    ax[1, 0].set_title("FV [ms-1]")
+    ax[1, 1].set_title("FT [m]")
+    # add colorbar
+    fig.colorbar(im3, ax=ax[1, 0])
+    fig.colorbar(im4, ax=ax[1, 1])
+
+    # add resistance area as transparent field
+    ax[1, 0].imshow(
+        np.flipud(np.where(fields["cResRasterOrig"] > 0, 1, np.nan)),
+        alpha=0.5,
+        extent=extentCellCenters,
+        zorder=4,
+    )
+    ax[1, 1].imshow(
+        np.flipud(np.where(fields["cResRasterOrig"] > 0, 1, np.nan)),
+        alpha=0.5,
+        extent=extentCellCenters,
+        zorder=4,
+    )
+
+    # fetch colormaps for FV and FT colorcode values below and above thresholds, mask 0 values
+    dmMin, dmMax = setMinMax(dmDet)
+    FTDetMin, FTDetMax = setMinMax(FTDet)
+
+    cmap5, col, ticks, norm = pU.makeColorMap(pU.colorMaps["dmDet"], dmMin, dmMax, continuous=pU.contCmap)
+    cmap6, col, ticks, norm = pU.makeColorMap(
+        pU.colorMaps["FTDet"], FTDetMin, FTDetMax, continuous=pU.contCmap
+    )
+    cmap5.set_bad(color="k", alpha=0.0)
+    cmap6.set_bad(color="k", alpha=0.0)
+    # add FV and FT fields
+    im5 = ax[0, 2].imshow(dmDet, cmap=cmap5, vmin=dmMin, vmax=dmMax, extent=extentCellCenters, zorder=4)
+    im6 = ax[1, 2].imshow(
+        FTDet,
+        cmap=cmap6,
+        vmin=FTDetMin,
+        vmax=FTDetMax,
+        extent=extentCellCenters,
+        zorder=4,
+    )
+    # set title
+    ax[0, 2].set_title("dmDet [kg]")
+    ax[1, 2].set_title("FTDet [m]")
+    # add colorbar
+    fig.colorbar(im5, ax=ax[0, 2])
+    fig.colorbar(im6, ax=ax[1, 2])
+
+    # save and or plot
+    plotName = "resistance_t%04d" % (tPlot)
+    outDir = pathlib.Path(cfg["avalancheDir"], "Outputs", "com1DFA", "cResPlots")
+    fU.makeADir(outDir)
+    pU.saveAndOrPlot({"pathResult": outDir}, plotName, fig)
+
+
+def setMinMax(resTypeField):
+    if np.isnan(resTypeField).all():
+        minVal = 0
+        maxVal = 0
+    else:
+        minVal = np.nanmin(resTypeField)
+        maxVal = np.nanmax(resTypeField)
+
+    return minVal, maxVal
+
+
+def massPlot(infoDict, massDetrainedTotal, t, avaDir, logName):
+    massDetrained = infoDict["massDetrained"]
+    massTotal = infoDict["massTotal"]
+    pfvTimeMax = infoDict["pfvTimeMax"]
+
+    # create figure
+    fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(pU.figW, pU.figH))
+
+    l1 = ax.plot(t, (100.0 / massTotal[0]) * np.asarray(massTotal), "k-", label="total mass")
+    ax2 = ax.twinx()
+    l2 = ax2.plot(
+        t,
+        (100.0 / massTotal[0]) * np.asarray(massDetrained),
+        color="royalblue",
+        linestyle="--",
+        label="detrained mass",
+    )
+    l3 = ax.plot(
+        t,
+        (100.0 / massTotal[0]) * (np.asarray(massTotal) - massDetrainedTotal),
+        "r--",
+        label="total mass + detrained mass",
+    )
+    ax3 = ax.twinx()
+    ax3.plot(t, pfvTimeMax, color="lightgray", linestyle="-")
+
+    ax.set_xlabel("time [s]")
+    ax.set_ylabel("total mass [% of initial mass]")
+    ax2.set_ylabel("detrained mass [% of initial mass]", color="royalblue")
+    ax3.set_ylabel("max velocity [ms-1]", color="lightgray")
+
+    ax2.spines["right"].set_color("royalblue")
+    ax2.tick_params(axis="y", colors="royalblue")
+    ax3.spines["right"].set_color("lightgray")
+    ax3.tick_params(axis="y", colors="lightgray")
+    ax3.spines["right"].set_position(("axes", 1.2))
+
+    lns = l1 + l2 + l3
+    labs = [l.get_label() for l in lns]
+    plt.legend(lns, labs, fontsize=9)
+
+    # save and or plot
+    plotName = "detrainedMass_%s" % (logName)
+    outDir = pathlib.Path(avaDir, "Outputs", "com1DFA", "reports")
     fU.makeADir(outDir)
     pU.saveAndOrPlot({"pathResult": outDir}, plotName, fig)
diff --git a/avaframe/out3Plot/plotUtils.py b/avaframe/out3Plot/plotUtils.py
index 03cbc7b65..740149d50 100644
--- a/avaframe/out3Plot/plotUtils.py
+++ b/avaframe/out3Plot/plotUtils.py
@@ -1,7 +1,7 @@
 """
-    Plot settings for output figures
+Plot settings for output figures
 
-    This file is part of Avaframe.
+This file is part of Avaframe.
 """
 
 import warnings
@@ -177,7 +177,7 @@
 # multi sequential colormap for Travel Angle
 levTA = list(fU.splitIniValueToArraySteps(cfgPlotUtils["travelAngleColorLevels"]))
 # Nuuk color map
-colorsTA = ['#fefeb2', '#d2d184', '#bab98d', '#a0a598', '#6f878d', '#386982', '#05598c']
+colorsTA = ["#fefeb2", "#d2d184", "#bab98d", "#a0a598", "#6f878d", "#386982", "#05598c"]
 cmapTA = copy.copy(cmapCrameri.nuuk.reversed())
 
 # colormap used if no resType provided
@@ -248,7 +248,8 @@
     "TA": cmapTravelAngle,
     "pke": cmapEnergy,
     "zdelta": cmapZdelta,
-    "dmDet": cmapProb
+    "dmDet": cmapProb,
+    "FTDet": cmapThickness,
 }
 
 cmapDEM = cmapGreys
@@ -396,8 +397,8 @@ def NonUnifIm(ax, x, y, z, xlab, ylab, aspect=None, **kwargs):
     ax.set_ylim([y.min(), y.max()])
     ax.set_xlabel(xlab)
     ax.set_ylabel(ylab)
-    if aspect == 'equal':
-        ax.set_aspect('equal')
+    if aspect == "equal":
+        ax.set_aspect("equal")
     return ref, im
 
 
@@ -910,7 +911,7 @@ def fetchContourCoords(xGrid, yGrid, data, level):
         # use codes of path
         labelID = 0
         if isinstance(contourP.get_paths()[0].codes, np.ndarray) == False:
-            log.warning('No contour found for level %.2f' % level)
+            log.warning("No contour found for level %.2f" % level)
         else:
             for ind, val in enumerate(contourP.get_paths()[0].codes):
                 if val == 1:
@@ -1037,7 +1038,14 @@ def createExtent(rowsMin, rowsMax, colsMin, colsMax, header, originLLCenter=True
     # create extent for cell centers lower left to upper right in meters
     extentCellCenters = [colsMinPlot, colsMaxPlot, rowsMinPlot, rowsMaxPlot]
 
-    return extentCellCenters, extentCellCorners, rowsMinPlot, rowsMaxPlot, colsMinPlot, colsMaxPlot
+    return (
+        extentCellCenters,
+        extentCellCorners,
+        rowsMinPlot,
+        rowsMaxPlot,
+        colsMinPlot,
+        colsMaxPlot,
+    )
 
 
 def createExtentMinMax(data, header, originLLCenter=True):
@@ -1064,22 +1072,22 @@ def createExtentMinMax(data, header, originLLCenter=True):
     nrows = data.shape[0]
     ncols = data.shape[1]
 
-    if nrows != header['nrows'] or ncols != header['ncols']:
-        message = 'Shape of data does not match nrows or ncols provided in header'
+    if nrows != header["nrows"] or ncols != header["ncols"]:
+        message = "Shape of data does not match nrows or ncols provided in header"
         log.error(message)
         raise AssertionError(message)
 
     # set extent in meters using cellSize and llcenter location
     if originLLCenter:
         rowsMinPlot = header["yllcenter"]
-        rowsMaxPlot = (nrows-1) * cellSize + header["yllcenter"]
+        rowsMaxPlot = (nrows - 1) * cellSize + header["yllcenter"]
         colsMinPlot = header["xllcenter"]
-        colsMaxPlot = (ncols-1) * cellSize + header["xllcenter"]
+        colsMaxPlot = (ncols - 1) * cellSize + header["xllcenter"]
     else:
         rowsMinPlot = 0
-        rowsMaxPlot = (nrows-1) * cellSize
+        rowsMaxPlot = (nrows - 1) * cellSize
         colsMinPlot = 0
-        colsMaxPlot = (ncols-1) * cellSize
+        colsMaxPlot = (ncols - 1) * cellSize
 
     # create extent of cell corners lower left to upper right in meters
     extentCellCorners = [
diff --git a/avaframe/out3Plot/plotUtilsCfg.ini b/avaframe/out3Plot/plotUtilsCfg.ini
index 4022723f4..ad0585e15 100644
--- a/avaframe/out3Plot/plotUtilsCfg.ini
+++ b/avaframe/out3Plot/plotUtilsCfg.ini
@@ -73,6 +73,7 @@ unitP = kPa
 unitFV = $ms^{-1}$
 unitrelth = m
 unitFT = m
+unitFTDet = m
 unitFM = kg
 unitFTV = $m^2s^{-1}$
 unitVx = $ms^{-1}$
@@ -115,6 +116,7 @@ contourLevelsFV = 0.5|1|5|10|25|50
 contourLevelsFT = 0.1|0.25|0.5|0.75|1
 contourLevelsTA = 10|15|20|25|30|35
 contourLevelsdmDet = -20|-15|-10|-5|0
+contourLevelsFTDet =  -0.1|-0.25|-0.5|-0.75|-1
 # name for result parameters
 nameppr = peak pressure
 namepft = peak flow thickness
@@ -124,6 +126,7 @@ namepke = peak kinetic energy
 nameRelTh = release thickness
 namedeltath = release thickness per particle
 nameFT = flow thickness
+nameFTDet = detrainment thickness
 nameFM = flow mass
 nameP = pressure
 nameFV = flow velocity
diff --git a/avaframe/tests/test_DFAfunctionsCython.py b/avaframe/tests/test_DFAfunctionsCython.py
index 1f01be1c8..210c7f2ec 100644
--- a/avaframe/tests/test_DFAfunctionsCython.py
+++ b/avaframe/tests/test_DFAfunctionsCython.py
@@ -1,4 +1,5 @@
-""
+""""""
+
 import numpy as np
 import math
 import configparser
@@ -11,75 +12,70 @@
 
 
 def test_getNeighborsC(capfd):
-    """ Test the grid search/particle location method"""
+    """Test the grid search/particle location method"""
     header = {}
-    header['ncols'] = 5
-    header['nrows'] = 6
-    header['cellsize'] = 1
+    header["ncols"] = 5
+    header["nrows"] = 6
+    header["cellsize"] = 1
     dem = {}
-    dem['header'] = header
-    dem['headerNeighbourGrid'] = header
+    dem["header"] = header
+    dem["headerNeighbourGrid"] = header
     particles = {}
-    particles['nPart'] = 18
-    particles['x'] = np.array(
-        [1.6, 0.4, 1, 2, 1, 2, 0, 1, 0, 2, 0, 2, 1, 2, 3, 3, 4, 0])
-    particles['y'] = np.array(
-        [2.6, 1.4, 0, 1, 3, 3, 2, 1, 0, 0, 3, 2, 2, 1, 1, 4, 5, 5])
-    particles['z'] = np.array(
-        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
-    particles['m'] = particles['z']
+    particles["nPart"] = 18
+    particles["x"] = np.array([1.6, 0.4, 1, 2, 1, 2, 0, 1, 0, 2, 0, 2, 1, 2, 3, 3, 4, 0])
+    particles["y"] = np.array([2.6, 1.4, 0, 1, 3, 3, 2, 1, 0, 0, 3, 2, 2, 1, 1, 4, 5, 5])
+    particles["z"] = np.array(
+        [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    )
+    particles["m"] = particles["z"]
     atol = 1e-10
-    indPCell = np.array([0.,  # always an extra zero at the begining
-                         1,  # found 1 particle
-                         2,  # found 1 particle
-                         3,  # found 1 particle
-                         3,  # nothing happens
-                         3,  # nothing happens
-                         4,  # found 1 particle
-                         5,  # found 1 particle
-                         7,   # found 2 particles
-                         8,  # found 1 particle
-                         8,  # nothing happens
-                         9,  # found 1 particle
-                         10,  # found 1 particle
-                         11,   # found 1 particles
-                         11, 11,  # nothing happens
-                         12,  # found 1 particle
-                         13,  # found 1 particle
-                         15,  # found 2 particle
-                         15, 15, 15, 15, 15,  # nothing happens
-                         16,  # found 1 particle
-                         16,  # nothing happens
-                         17,  # found 1 particle
-                         17, 17, 17,  # nothing happens
-                         18])  # found 1 particle
-    pInC = np.array([8,
-                     2,
-                     9,
-                     1,
-                     7,
-                     13, 3,
-                     14,
-                     6,
-                     12,
-                     11,
-                     10,
-                     4,
-                     5, 0,
-                     15,
-                     17,
-                     16])
+    indPCell = np.array(
+        [
+            0.0,  # always an extra zero at the begining
+            1,  # found 1 particle
+            2,  # found 1 particle
+            3,  # found 1 particle
+            3,  # nothing happens
+            3,  # nothing happens
+            4,  # found 1 particle
+            5,  # found 1 particle
+            7,  # found 2 particles
+            8,  # found 1 particle
+            8,  # nothing happens
+            9,  # found 1 particle
+            10,  # found 1 particle
+            11,  # found 1 particles
+            11,
+            11,  # nothing happens
+            12,  # found 1 particle
+            13,  # found 1 particle
+            15,  # found 2 particle
+            15,
+            15,
+            15,
+            15,
+            15,  # nothing happens
+            16,  # found 1 particle
+            16,  # nothing happens
+            17,  # found 1 particle
+            17,
+            17,
+            17,  # nothing happens
+            18,
+        ]
+    )  # found 1 particle
+    pInC = np.array([8, 2, 9, 1, 7, 13, 3, 14, 6, 12, 11, 10, 4, 5, 0, 15, 17, 16])
 
     particles = DFAfunC.getNeighborsC(particles, dem)
-#    print(particles['inCellDEM'])
-#    print(particles['indPartInCell'])
-#    print(particles['partInCell'])
-    assert np.allclose(particles['indPartInCell'], indPCell, atol=atol)
-    assert np.allclose(particles['partInCell'], pInC, atol=atol)
+    #    print(particles['inCellDEM'])
+    #    print(particles['indPartInCell'])
+    #    print(particles['partInCell'])
+    assert np.allclose(particles["indPartInCell"], indPCell, atol=atol)
+    assert np.allclose(particles["partInCell"], pInC, atol=atol)
 
 
 def test_computeEntMassAndForce(capfd):
-    """ Test the computeEntMassAndForce function"""
+    """Test the computeEntMassAndForce function"""
     dt = 0.1
     entrMassCell = 0
     areaPart = 1
@@ -88,29 +84,32 @@ def test_computeEntMassAndForce(capfd):
     entEroEnergy = 5000
     rhoEnt = 100
     dm, areaEntrPart = DFAfunC.computeEntMassAndForce(
-        dt, entrMassCell, areaPart, uMag, tau, entEroEnergy, rhoEnt)
-#    print(dm, areaEntrPart)
+        dt, entrMassCell, areaPart, uMag, tau, entEroEnergy, rhoEnt
+    )
+    #    print(dm, areaEntrPart)
     assert dm == 0
     assert areaEntrPart == 1
 
     entrMassCell = 200
     entEroEnergy = 0
     dm, areaEntrPart = DFAfunC.computeEntMassAndForce(
-        dt, entrMassCell, areaPart, uMag, tau, entEroEnergy, rhoEnt)
-#    print(dm, areaEntrPart)
+        dt, entrMassCell, areaPart, uMag, tau, entEroEnergy, rhoEnt
+    )
+    #    print(dm, areaEntrPart)
     assert dm == 200
     assert areaEntrPart == 2
 
     entEroEnergy = 5000
     dm, areaEntrPart = DFAfunC.computeEntMassAndForce(
-        dt, entrMassCell, areaPart, uMag, tau, entEroEnergy, rhoEnt)
-#    print(dm, areaEntrPart)
+        dt, entrMassCell, areaPart, uMag, tau, entEroEnergy, rhoEnt
+    )
+    #    print(dm, areaEntrPart)
     assert dm == 0.2
     assert areaEntrPart == 1
-    
+
 
 def test_computeDetMass(capfd):
-    """ Test the computeDetMass function"""
+    """Test the computeDetMass function"""
     dt = 0.1
     detCell = 0
     areaPart = 1
@@ -120,55 +119,32 @@ def test_computeDetMass(capfd):
 
     detCell = 10
     dmDet = DFAfunC.computeDetMass(dt, detCell, areaPart, uMag)
-#    print(dmDet)
+    #    print(dmDet)
     assert dmDet == -0.1
 
 
 def test_computeResForce(capfd):
-    """ Test the computeResForce function"""
-    hRes = 2
-    h = 1
+    """Test the computeResForce function"""
+
     areaPart = 1
     rho = 200
     cResCell = 1
     uMag = 10
     explicitFriction = 0
-    #cRes
+    # cRes
     resistanceType = 1
-    cResPart = DFAfunC.computeResForce(
-        hRes, h, areaPart, rho, cResCell, uMag, explicitFriction, resistanceType)
+    cResPart = DFAfunC.computeResForce(areaPart, rho, cResCell, uMag, explicitFriction, resistanceType)
     print(cResPart)
     assert cResPart == -2000
 
-    h = 3
-    cResPart = DFAfunC.computeResForce(
-        hRes, h, areaPart, rho, cResCell, uMag, explicitFriction, resistanceType)
-    print(cResPart)
-    assert cResPart == -4000
-
     explicitFriction = 1
-    cResPart = DFAfunC.computeResForce(
-        hRes, h, areaPart, rho, cResCell, uMag, explicitFriction, resistanceType)
+    cResPart = DFAfunC.computeResForce(areaPart, rho, cResCell, uMag, explicitFriction, resistanceType)
     print(cResPart)
-    assert cResPart == -40000
-
-    #cResH
-    explicitFriction = 0
-    resistanceType = 2
-    cResPart = DFAfunC.computeResForce(
-        hRes, h, areaPart, rho, cResCell, uMag, explicitFriction, resistanceType)
-    print(cResPart)
-    assert cResPart == -2000
-
-    h = 1
-    cResPart = DFAfunC.computeResForce(
-        hRes, h, areaPart, rho, cResCell, uMag, explicitFriction, resistanceType)
-    print(cResPart)
-    assert cResPart == -2000
+    assert cResPart == -20000
 
 
 def test_account4FrictionForce(capfd):
-    """ Test the account4FrictionForce function"""
+    """Test the account4FrictionForce function"""
     uxNew = 10
     uyNew = 0
     uzNew = 0
@@ -178,8 +154,9 @@ def test_account4FrictionForce(capfd):
     uMag = 10
     explicitFriction = 0
     uxNew, uyNew, uzNew, dtStop = DFAfunC.account4FrictionForce(
-        uxNew, uyNew, uzNew, m, dt, forceFrict, uMag, explicitFriction)
-#    print(uxNew, uyNew, uzNew, dtStop)
+        uxNew, uyNew, uzNew, m, dt, forceFrict, uMag, explicitFriction
+    )
+    #    print(uxNew, uyNew, uzNew, dtStop)
     assert dtStop == 0.1
     assert uxNew == 5
 
@@ -189,9 +166,10 @@ def test_account4FrictionForce(capfd):
     explicitFriction = 1
     m = 0.5
     uxNew, uyNew, uzNew, dtStop = DFAfunC.account4FrictionForce(
-        uxNew, uyNew, uzNew, m, dt, forceFrict, uMag, explicitFriction)
-#    print(uxNew, uyNew, uzNew, dtStop)
-#    print(dt*forceFrict/m)
+        uxNew, uyNew, uzNew, m, dt, forceFrict, uMag, explicitFriction
+    )
+    #    print(uxNew, uyNew, uzNew, dtStop)
+    #    print(dt*forceFrict/m)
     assert dtStop == 0.05
     assert uxNew == 0
 
@@ -200,192 +178,307 @@ def test_account4FrictionForce(capfd):
     uzNew = 0
     m = 10
     uxNew, uyNew, uzNew, dtStop = DFAfunC.account4FrictionForce(
-        uxNew, uyNew, uzNew, m, dt, forceFrict, uMag, explicitFriction)
-#    print(uxNew, uyNew, uzNew, dtStop)
-#    print(dt*forceFrict/m)
+        uxNew, uyNew, uzNew, m, dt, forceFrict, uMag, explicitFriction
+    )
+    #    print(uxNew, uyNew, uzNew, dtStop)
+    #    print(dt*forceFrict/m)
     assert dtStop == 0.1
     assert uxNew == 9.0
 
 
 def test_updatePositionC():
-    """ test updating the position of particles """
+    """test updating the position of particles"""
 
     # TODO: make test also if velocity in z not zero!! - so to test also when reprojecting onto surface
 
     # setup required input
     cfg = configparser.ConfigParser()
-    cfg['GENERAL'] = {'stopCrit': '0.01', 'stopCritIni': '0.1', 'stopCritIniSmall': '1.001', 'stopCritType': 'kinEnergy',
-                      'uFlowingThreshold': '0.1', 'gravAcc': '9.81', 'velMagMin': '1.e-6',  'rho': '100.',
-                      'interpOption': '2',   'explicitFriction': '0', 'centeredPosition': '1',
-                      'reprojMethodPosition': '2', 'reprojectionIterations': '5', 'thresholdProjection': '0.001',
-                      'dissDam': '1', 'snowSlide': '1', 'wetSnow': '1'}
-
-    particles = {'dt': 1.0, 'm': np.asarray([10., 10., 10.]), 'idFixed': np.asarray([0., 0., 0.]), 'trajectoryLengthXY': np.asarray([0., 0., 0.]),
-                  'trajectoryLengthXYCor': np.asarray([0., 0., 0.]), 'trajectoryLengthXYZ': np.asarray([0., 0., 0.]), 'x': np.asarray([0., 1., 2.]), 'y': np.asarray([2., 3., 4.]),
-                  'z': np.asarray([1., 1., 1.]), 'ux': np.asarray([1., 1., 1.]), 'uy': np.asarray([1., 1., 1.]),
-                  'uz': np.asarray([0., 0., 0.]), 'uAcc': np.asarray([0., 0., 0.]), 'kineticEne': 0.0, 'peakKinEne': 0.0,
-                  'peakForceSPH': 0.0, 'forceSPHIni': 0.0, 'nPart': 3,
-                  'peakMassFlowing': 0.0, 'iterate': True, 'totalEnthalpy': np.asarray([0., 0., 0.]),
-                  'velocityMag': np.asarray([1., 1., 1.])}
-    particles['potentialEne'] = np.sum(9.81 * particles['z'] * particles['m'])
+    cfg["GENERAL"] = {
+        "stopCrit": "0.01",
+        "stopCritIni": "0.1",
+        "stopCritIniSmall": "1.001",
+        "stopCritType": "kinEnergy",
+        "uFlowingThreshold": "0.1",
+        "gravAcc": "9.81",
+        "velMagMin": "1.e-6",
+        "rho": "100.",
+        "interpOption": "2",
+        "explicitFriction": "0",
+        "centeredPosition": "1",
+        "reprojMethodPosition": "2",
+        "reprojectionIterations": "5",
+        "thresholdProjection": "0.001",
+        "dissDam": "1",
+        "snowSlide": "1",
+        "wetSnow": "1",
+    }
+
+    particles = {
+        "dt": 1.0,
+        "m": np.asarray([10.0, 10.0, 10.0]),
+        "idFixed": np.asarray([0.0, 0.0, 0.0]),
+        "trajectoryLengthXY": np.asarray([0.0, 0.0, 0.0]),
+        "trajectoryLengthXYCor": np.asarray([0.0, 0.0, 0.0]),
+        "trajectoryLengthXYZ": np.asarray([0.0, 0.0, 0.0]),
+        "x": np.asarray([0.0, 1.0, 2.0]),
+        "y": np.asarray([2.0, 3.0, 4.0]),
+        "z": np.asarray([1.0, 1.0, 1.0]),
+        "ux": np.asarray([1.0, 1.0, 1.0]),
+        "uy": np.asarray([1.0, 1.0, 1.0]),
+        "uz": np.asarray([0.0, 0.0, 0.0]),
+        "uAcc": np.asarray([0.0, 0.0, 0.0]),
+        "kineticEne": 0.0,
+        "peakKinEne": 0.0,
+        "peakForceSPH": 0.0,
+        "forceSPHIni": 0.0,
+        "nPart": 3,
+        "peakMassFlowing": 0.0,
+        "iterate": True,
+        "totalEnthalpy": np.asarray([0.0, 0.0, 0.0]),
+        "velocityMag": np.asarray([1.0, 1.0, 1.0]),
+    }
+    particles["potentialEne"] = np.sum(9.81 * particles["z"] * particles["m"])
 
     demHeader = {}
-    demHeader['xllcenter'] = 0.0
-    demHeader['yllcenter'] = 0.0
-    demHeader['cellsize'] = 5.0
-    demHeader['nodata_value'] = -9999
-    demHeader['nrows'] = 10
-    demHeader['ncols'] = 10
-    dem = {'header': demHeader}
-    dem['rasterData'] = np.ones((demHeader['nrows'], demHeader['ncols']))
-    dem['outOfDEM'] = np.zeros((demHeader['nrows'], demHeader['ncols'])).flatten()
-    dem['Nx'] = np.zeros((demHeader['nrows'], demHeader['ncols']))
-    dem['Ny'] = np.zeros((demHeader['nrows'], demHeader['ncols']))
-    dem['Nz'] = np.ones((demHeader['nrows'], demHeader['ncols']))
-
-    force = {'forceZ': np.asarray([0., 0., 0.]), 'forceFrict': np.asarray([10., 10., 10.]),
-             'dM': np.asarray([0., 0., 0.]), 'dMDet': np.asarray([0., 0., 0.]),
-             'forceX': np.asarray([50., 50., 50.]), 'forceY': np.asarray([50., 50., 50.]),
-             'forceSPHX': np.asarray([50., 50., 50.]),
-             'forceSPHY': np.asarray([50., 50., 50.]), 'forceSPHZ': np.asarray([0., 0., 0.])}
-    fields = {'FT': np.zeros((2, 2))}
+    demHeader["xllcenter"] = 0.0
+    demHeader["yllcenter"] = 0.0
+    demHeader["cellsize"] = 5.0
+    demHeader["nodata_value"] = -9999
+    demHeader["nrows"] = 10
+    demHeader["ncols"] = 10
+    dem = {"header": demHeader}
+    dem["rasterData"] = np.ones((demHeader["nrows"], demHeader["ncols"]))
+    dem["outOfDEM"] = np.zeros((demHeader["nrows"], demHeader["ncols"])).flatten()
+    dem["Nx"] = np.zeros((demHeader["nrows"], demHeader["ncols"]))
+    dem["Ny"] = np.zeros((demHeader["nrows"], demHeader["ncols"]))
+    dem["Nz"] = np.ones((demHeader["nrows"], demHeader["ncols"]))
+
+    force = {
+        "forceZ": np.asarray([0.0, 0.0, 0.0]),
+        "forceFrict": np.asarray([10.0, 10.0, 10.0]),
+        "dM": np.asarray([0.0, 0.0, 0.0]),
+        "dMDet": np.asarray([0.0, 0.0, 0.0]),
+        "forceX": np.asarray([50.0, 50.0, 50.0]),
+        "forceY": np.asarray([50.0, 50.0, 50.0]),
+        "forceSPHX": np.asarray([50.0, 50.0, 50.0]),
+        "forceSPHY": np.asarray([50.0, 50.0, 50.0]),
+        "forceSPHZ": np.asarray([0.0, 0.0, 0.0]),
+    }
+    fields = {"FT": np.zeros((2, 2))}
     # crete a dummy dict (needed so that cython runs)
-    wallLineDict = {'dam': 0, 'nIterDam': 1, 'cellsCrossed': np.zeros((dem['header']['ncols']*dem['header']['nrows'])).astype(int)}
-    for key in ['x', 'y', 'z', 'xCrown', 'yCrown', 'zCrown', 'xTangent', 'yTangent', 'zTangent']:
-        wallLineDict[key] = np.ones((1))*1.0
-    for key in ['nPoints', 'height', 'slope', 'restitutionCoefficient']:
+    wallLineDict = {
+        "dam": 0,
+        "nIterDam": 1,
+        "cellsCrossed": np.zeros((dem["header"]["ncols"] * dem["header"]["nrows"])).astype(int),
+    }
+    for key in [
+        "x",
+        "y",
+        "z",
+        "xCrown",
+        "yCrown",
+        "zCrown",
+        "xTangent",
+        "yTangent",
+        "zTangent",
+    ]:
+        wallLineDict[key] = np.ones((1)) * 1.0
+    for key in ["nPoints", "height", "slope", "restitutionCoefficient"]:
         wallLineDict[key] = 0
-    dem['damLine'] = wallLineDict
+    dem["damLine"] = wallLineDict
     typeStop = 0
 
     # kinetic energy new
     kinEneNew = 0.0
     potEneNew = 0.0
     for k in range(3):
-        kinEneNew = kinEneNew + particles['m'][k] * np.sqrt(5.5**2 + 5.5**2 + 0**2)**2 * 0.5
-        potEneNew = potEneNew + particles['m'][k] * 9.81 + 0.0
-
-    particles = DFAfunC.updatePositionC(cfg['GENERAL'], particles, dem, force, fields, typeStop=typeStop)
-    uAcc = (np.sqrt(5.5**2 + 5.5**2 +0.) - np.sqrt(1.**2 + 1.**2 + 0.)) / 1.
-    velocityMag = np.asarray([np.sqrt((5.5**2) + (5.5**2)), np.sqrt((5.5**2) + (5.5**2)), np.sqrt((5.5**2) + (5.5**2))])
-
-    assert np.array_equal(particles['m'], np.asarray([10., 10., 10.]))
-    assert np.array_equal(particles['x'], np.array([3.25, 4.25, 5.25]))
-    assert np.array_equal(particles['y'], np.asarray([5.25, 6.25, 7.25]))
-    assert np.array_equal(particles['z'], np.asarray([1., 1., 1.]))
-    assert np.allclose(particles['ux'], np.asarray([5.5, 5.5, 5.5]), atol=1.e-4)
-    assert np.allclose(particles['uy'], np.asarray([5.5, 5.5, 5.5]), atol=1.e-4)
-    assert np.allclose(particles['uAcc'], np.asarray([uAcc, uAcc, uAcc]), atol=1.e-4)
-    assert np.allclose(particles['velocityMag'], velocityMag, atol=1.e-4)
-    assert np.array_equal(particles['uz'], np.asarray([0., 0., 0.]))
-    assert particles['massEntrained'] == 0.0
-    assert particles['massDetrained'] == 0.0
-    assert particles['nPart'] == 3
-    assert (kinEneNew- 1.e-4) < particles['kineticEne'] < (kinEneNew+1.e-4)
-    assert (potEneNew-1.e-4) < particles['potentialEne'] < (potEneNew +1.e-4)
-    assert particles['iterate'] == True
-
-    particles = {'dt': 1.0, 'm': np.asarray([10., 10., 10.]), 'idFixed': np.asarray([0., 0., 0.]), 'trajectoryLengthXY': np.asarray([0., 0., 0.]),
-                  'trajectoryLengthXYCor': np.asarray([0., 0., 0.]), 'trajectoryLengthXYZ': np.asarray([0., 0., 0.]), 'x': np.asarray([0., 1., 2.]), 'y': np.asarray([2., 3., 4.]),
-                  'z': np.asarray([1., 1., 1.]), 'ux': np.asarray([1., 1., 1.]), 'uy': np.asarray([1., 1., 1.]),
-                  'uz': np.asarray([0., 0., 0.]), 'uAcc': np.asarray([0., 0., 0.]), 'kineticEne': 0.0, 'peakKinEne': 100000.0,
-                  'peakForceSPH': 0.0, 'forceSPHIni': 0.0, 'nPart': 3,
-                  'velocityMag': np.asarray([1., 1., 1.]),
-                  'peakMassFlowing': 0.0, 'iterate': True, 'totalEnthalpy': np.asarray([0., 0., 0.])}
-    particles['potentialEne'] = np.sum(9.81 * particles['z'] * particles['m'])
+        kinEneNew = kinEneNew + particles["m"][k] * np.sqrt(5.5**2 + 5.5**2 + 0**2) ** 2 * 0.5
+        potEneNew = potEneNew + particles["m"][k] * 9.81 + 0.0
+
+    particles = DFAfunC.updatePositionC(cfg["GENERAL"], particles, dem, force, fields, typeStop=typeStop)
+    uAcc = (np.sqrt(5.5**2 + 5.5**2 + 0.0) - np.sqrt(1.0**2 + 1.0**2 + 0.0)) / 1.0
+    velocityMag = np.asarray(
+        [
+            np.sqrt((5.5**2) + (5.5**2)),
+            np.sqrt((5.5**2) + (5.5**2)),
+            np.sqrt((5.5**2) + (5.5**2)),
+        ]
+    )
+
+    assert np.array_equal(particles["m"], np.asarray([10.0, 10.0, 10.0]))
+    assert np.array_equal(particles["x"], np.array([3.25, 4.25, 5.25]))
+    assert np.array_equal(particles["y"], np.asarray([5.25, 6.25, 7.25]))
+    assert np.array_equal(particles["z"], np.asarray([1.0, 1.0, 1.0]))
+    assert np.allclose(particles["ux"], np.asarray([5.5, 5.5, 5.5]), atol=1.0e-4)
+    assert np.allclose(particles["uy"], np.asarray([5.5, 5.5, 5.5]), atol=1.0e-4)
+    assert np.allclose(particles["uAcc"], np.asarray([uAcc, uAcc, uAcc]), atol=1.0e-4)
+    assert np.allclose(particles["velocityMag"], velocityMag, atol=1.0e-4)
+    assert np.array_equal(particles["uz"], np.asarray([0.0, 0.0, 0.0]))
+    assert particles["massEntrained"] == 0.0
+    assert particles["massDetrained"] == 0.0
+    assert particles["nPart"] == 3
+    assert (kinEneNew - 1.0e-4) < particles["kineticEne"] < (kinEneNew + 1.0e-4)
+    assert (potEneNew - 1.0e-4) < particles["potentialEne"] < (potEneNew + 1.0e-4)
+    assert particles["iterate"] == True
+
+    particles = {
+        "dt": 1.0,
+        "m": np.asarray([10.0, 10.0, 10.0]),
+        "idFixed": np.asarray([0.0, 0.0, 0.0]),
+        "trajectoryLengthXY": np.asarray([0.0, 0.0, 0.0]),
+        "trajectoryLengthXYCor": np.asarray([0.0, 0.0, 0.0]),
+        "trajectoryLengthXYZ": np.asarray([0.0, 0.0, 0.0]),
+        "x": np.asarray([0.0, 1.0, 2.0]),
+        "y": np.asarray([2.0, 3.0, 4.0]),
+        "z": np.asarray([1.0, 1.0, 1.0]),
+        "ux": np.asarray([1.0, 1.0, 1.0]),
+        "uy": np.asarray([1.0, 1.0, 1.0]),
+        "uz": np.asarray([0.0, 0.0, 0.0]),
+        "uAcc": np.asarray([0.0, 0.0, 0.0]),
+        "kineticEne": 0.0,
+        "peakKinEne": 100000.0,
+        "peakForceSPH": 0.0,
+        "forceSPHIni": 0.0,
+        "nPart": 3,
+        "velocityMag": np.asarray([1.0, 1.0, 1.0]),
+        "peakMassFlowing": 0.0,
+        "iterate": True,
+        "totalEnthalpy": np.asarray([0.0, 0.0, 0.0]),
+    }
+    particles["potentialEne"] = np.sum(9.81 * particles["z"] * particles["m"])
 
     # call function to be tested
-    particles = DFAfunC.updatePositionC(cfg['GENERAL'], particles, dem, force, fields, typeStop=typeStop)
-
-    assert np.array_equal(particles['m'], np.asarray([10., 10., 10.]))
-    assert np.array_equal(particles['x'], np.array([3.25, 4.25, 5.25]))
-    assert np.array_equal(particles['y'], np.asarray([5.25, 6.25, 7.25]))
-    assert np.array_equal(particles['z'], np.asarray([1., 1., 1.]))
-    assert np.allclose(particles['ux'], np.asarray([5.5, 5.5, 5.5]), atol=1.e-4)
-    assert np.allclose(particles['uy'], np.asarray([5.5, 5.5, 5.5]), atol=1.e-4)
-    assert np.allclose(particles['uAcc'], np.asarray([uAcc, uAcc, uAcc]), atol=1.e-4)
-    assert np.array_equal(particles['uz'], np.asarray([0., 0., 0.]))
-    assert particles['massEntrained'] == 0.0
-    assert particles['nPart'] == 3
-    assert (kinEneNew- 1.e-4) < particles['kineticEne'] < (kinEneNew+1.e-4)
-    assert (potEneNew-1.e-4) < particles['potentialEne'] < (potEneNew +1.e-4)
-    assert particles['iterate'] == False
-
-    particles = {'dt': 1.0, 'm': np.asarray([10., 10., 10.]), 'idFixed': np.asarray([0., 0., 0.]), 'trajectoryLengthXY': np.asarray([0., 0., 0.]),
-                  'trajectoryLengthXYCor': np.asarray([0., 0., 0.]), 'trajectoryLengthXYZ': np.asarray([0., 0., 0.]), 'x': np.asarray([0., 1., 2.]), 'y': np.asarray([2., 3., 4.]),
-                  'z': np.asarray([1., 1., 1.]), 'ux': np.asarray([1., 1., 1.]), 'uy': np.asarray([1., 1., 1.]),
-                  'uz': np.asarray([0., 0., 0.]), 'uAcc': np.asarray([0., 0., 0.]), 'kineticEne': 0.0, 'peakKinEne': 10000.0,
-                  'peakForceSPH': 100000.0, 'forceSPHIni': 1.e5, 'nPart': 3,
-                  'velocityMag': np.asarray([1., 1., 1.]),
-                  'peakMassFlowing': 0.0, 'iterate': True,  'totalEnthalpy': np.asarray([0., 0., 0.])}
-    particles['potentialEne'] = np.sum(9.81 * particles['z'] * particles['m'])
+    particles = DFAfunC.updatePositionC(cfg["GENERAL"], particles, dem, force, fields, typeStop=typeStop)
+
+    assert np.array_equal(particles["m"], np.asarray([10.0, 10.0, 10.0]))
+    assert np.array_equal(particles["x"], np.array([3.25, 4.25, 5.25]))
+    assert np.array_equal(particles["y"], np.asarray([5.25, 6.25, 7.25]))
+    assert np.array_equal(particles["z"], np.asarray([1.0, 1.0, 1.0]))
+    assert np.allclose(particles["ux"], np.asarray([5.5, 5.5, 5.5]), atol=1.0e-4)
+    assert np.allclose(particles["uy"], np.asarray([5.5, 5.5, 5.5]), atol=1.0e-4)
+    assert np.allclose(particles["uAcc"], np.asarray([uAcc, uAcc, uAcc]), atol=1.0e-4)
+    assert np.array_equal(particles["uz"], np.asarray([0.0, 0.0, 0.0]))
+    assert particles["massEntrained"] == 0.0
+    assert particles["nPart"] == 3
+    assert (kinEneNew - 1.0e-4) < particles["kineticEne"] < (kinEneNew + 1.0e-4)
+    assert (potEneNew - 1.0e-4) < particles["potentialEne"] < (potEneNew + 1.0e-4)
+    assert particles["iterate"] == False
+
+    particles = {
+        "dt": 1.0,
+        "m": np.asarray([10.0, 10.0, 10.0]),
+        "idFixed": np.asarray([0.0, 0.0, 0.0]),
+        "trajectoryLengthXY": np.asarray([0.0, 0.0, 0.0]),
+        "trajectoryLengthXYCor": np.asarray([0.0, 0.0, 0.0]),
+        "trajectoryLengthXYZ": np.asarray([0.0, 0.0, 0.0]),
+        "x": np.asarray([0.0, 1.0, 2.0]),
+        "y": np.asarray([2.0, 3.0, 4.0]),
+        "z": np.asarray([1.0, 1.0, 1.0]),
+        "ux": np.asarray([1.0, 1.0, 1.0]),
+        "uy": np.asarray([1.0, 1.0, 1.0]),
+        "uz": np.asarray([0.0, 0.0, 0.0]),
+        "uAcc": np.asarray([0.0, 0.0, 0.0]),
+        "kineticEne": 0.0,
+        "peakKinEne": 10000.0,
+        "peakForceSPH": 100000.0,
+        "forceSPHIni": 1.0e5,
+        "nPart": 3,
+        "velocityMag": np.asarray([1.0, 1.0, 1.0]),
+        "peakMassFlowing": 0.0,
+        "iterate": True,
+        "totalEnthalpy": np.asarray([0.0, 0.0, 0.0]),
+    }
+    particles["potentialEne"] = np.sum(9.81 * particles["z"] * particles["m"])
     typeStop = 1
 
     sphForceNew = 0.0
     kinEneNew = 0.0
     potEneNew = 0.0
     for k in range(3):
-        sphForceNew = sphForceNew + particles['m'][k] * np.sqrt(50.**2 +50.**2 + 0.**2)**2.
-        kinEneNew = kinEneNew + particles['m'][k] * np.sqrt(11.**2 +11.**2 + 0**2)**2 * 0.5
-        potEneNew = potEneNew + particles['m'][k] * 9.81 + 0.0
+        sphForceNew = sphForceNew + particles["m"][k] * np.sqrt(50.0**2 + 50.0**2 + 0.0**2) ** 2.0
+        kinEneNew = kinEneNew + particles["m"][k] * np.sqrt(11.0**2 + 11.0**2 + 0**2) ** 2 * 0.5
+        potEneNew = potEneNew + particles["m"][k] * 9.81 + 0.0
 
     # call function to be tested
-    particles = DFAfunC.updatePositionC(cfg['GENERAL'], particles, dem, force, fields, typeStop=typeStop)
-#    print('sph', particles['peakForceSPH'], sphForceNew)
-
-    uAcc = (np.sqrt(11.**2 + 11.0**2 +0.) - np.sqrt(1.**2 + 1.**2 + 0.)) / 1.
-    velocityMag = np.asarray([np.sqrt((11.**2) + (11.**2)), np.sqrt((11.**2) + (11.**2)), np.sqrt((11.**2) + (11.**2))])
-
-    assert np.array_equal(particles['m'], np.asarray([10., 10., 10.]))
-    assert np.array_equal(particles['x'], np.array([6., 7., 8.]))
-    assert np.array_equal(particles['y'], np.asarray([8., 9., 10.]))
-    assert np.array_equal(particles['z'], np.asarray([1., 1., 1.]))
-    assert np.allclose(particles['ux'], np.asarray([11., 11., 11.]), atol=1.e-4)
-    assert np.allclose(particles['uy'], np.asarray([11., 11., 11.]), atol=1.e-4)
-    assert np.allclose(particles['velocityMag'], velocityMag, atol=1.e-4)
-    assert np.array_equal(particles['uz'], np.asarray([0., 0., 0.]))
-    assert np.allclose(particles['uAcc'], np.asarray([uAcc, uAcc, uAcc]), atol=1.e-4)
-    assert particles['massEntrained'] == 0.0
-    assert particles['nPart'] == 3
-    assert (kinEneNew- 1.e-4) < particles['kineticEne'] < (kinEneNew+1.e-4)
-    assert (potEneNew-1.e-4) < particles['potentialEne'] < (potEneNew +1.e-4)
-    assert particles['iterate'] == False
-
-    particles = {'dt': 1.0, 'm': np.asarray([10., 10., 10.]), 'idFixed': np.asarray([0., 0., 0.]), 'trajectoryLengthXY': np.asarray([0., 0., 0.]),
-                  'trajectoryLengthXYCor': np.asarray([0., 0., 0.]), 'trajectoryLengthXYZ': np.asarray([0., 0., 0.]), 'x': np.asarray([0., 1., 2.]), 'y': np.asarray([2., 3., 4.]),
-                  'z': np.asarray([1., 1., 1.]), 'ux': np.asarray([1., 1., 1.]), 'uy': np.asarray([1., 1., 1.]),
-                  'uz': np.asarray([0., 0., 0.]), 'uAcc': np.asarray([0., 0., 0.]), 'kineticEne': 0.0, 'peakKinEne': 10000.0,
-                  'peakForceSPH': 1000.0, 'forceSPHIni': 1.e5, 'nPart': 3,
-                  'velocityMag': np.asarray([1., 1., 1.]),
-                  'peakMassFlowing': 0.0, 'iterate': True, 'totalEnthalpy': np.asarray([0., 0., 0.])}
-    particles['potentialEne'] = np.sum(9.81 * particles['z'] * particles['m'])
+    particles = DFAfunC.updatePositionC(cfg["GENERAL"], particles, dem, force, fields, typeStop=typeStop)
+    #    print('sph', particles['peakForceSPH'], sphForceNew)
+
+    uAcc = (np.sqrt(11.0**2 + 11.0**2 + 0.0) - np.sqrt(1.0**2 + 1.0**2 + 0.0)) / 1.0
+    velocityMag = np.asarray(
+        [
+            np.sqrt((11.0**2) + (11.0**2)),
+            np.sqrt((11.0**2) + (11.0**2)),
+            np.sqrt((11.0**2) + (11.0**2)),
+        ]
+    )
+
+    assert np.array_equal(particles["m"], np.asarray([10.0, 10.0, 10.0]))
+    assert np.array_equal(particles["x"], np.array([6.0, 7.0, 8.0]))
+    assert np.array_equal(particles["y"], np.asarray([8.0, 9.0, 10.0]))
+    assert np.array_equal(particles["z"], np.asarray([1.0, 1.0, 1.0]))
+    assert np.allclose(particles["ux"], np.asarray([11.0, 11.0, 11.0]), atol=1.0e-4)
+    assert np.allclose(particles["uy"], np.asarray([11.0, 11.0, 11.0]), atol=1.0e-4)
+    assert np.allclose(particles["velocityMag"], velocityMag, atol=1.0e-4)
+    assert np.array_equal(particles["uz"], np.asarray([0.0, 0.0, 0.0]))
+    assert np.allclose(particles["uAcc"], np.asarray([uAcc, uAcc, uAcc]), atol=1.0e-4)
+    assert particles["massEntrained"] == 0.0
+    assert particles["nPart"] == 3
+    assert (kinEneNew - 1.0e-4) < particles["kineticEne"] < (kinEneNew + 1.0e-4)
+    assert (potEneNew - 1.0e-4) < particles["potentialEne"] < (potEneNew + 1.0e-4)
+    assert particles["iterate"] == False
+
+    particles = {
+        "dt": 1.0,
+        "m": np.asarray([10.0, 10.0, 10.0]),
+        "idFixed": np.asarray([0.0, 0.0, 0.0]),
+        "trajectoryLengthXY": np.asarray([0.0, 0.0, 0.0]),
+        "trajectoryLengthXYCor": np.asarray([0.0, 0.0, 0.0]),
+        "trajectoryLengthXYZ": np.asarray([0.0, 0.0, 0.0]),
+        "x": np.asarray([0.0, 1.0, 2.0]),
+        "y": np.asarray([2.0, 3.0, 4.0]),
+        "z": np.asarray([1.0, 1.0, 1.0]),
+        "ux": np.asarray([1.0, 1.0, 1.0]),
+        "uy": np.asarray([1.0, 1.0, 1.0]),
+        "uz": np.asarray([0.0, 0.0, 0.0]),
+        "uAcc": np.asarray([0.0, 0.0, 0.0]),
+        "kineticEne": 0.0,
+        "peakKinEne": 10000.0,
+        "peakForceSPH": 1000.0,
+        "forceSPHIni": 1.0e5,
+        "nPart": 3,
+        "velocityMag": np.asarray([1.0, 1.0, 1.0]),
+        "peakMassFlowing": 0.0,
+        "iterate": True,
+        "totalEnthalpy": np.asarray([0.0, 0.0, 0.0]),
+    }
+    particles["potentialEne"] = np.sum(9.81 * particles["z"] * particles["m"])
     typeStop = 1
 
     sphForceNew = 0.0
     kinEneNew = 0.0
     potEneNew = 0.0
     for k in range(3):
-        sphForceNew = sphForceNew + particles['m'][k] * np.sqrt(50.**2 +50.**2 + 0.**2)**2.
-        kinEneNew = kinEneNew + particles['m'][k] * np.sqrt(11.**2 +11.**2 + 0**2)**2 * 0.5
-        potEneNew = potEneNew + particles['m'][k] * 9.81 + 0.0
+        sphForceNew = sphForceNew + particles["m"][k] * np.sqrt(50.0**2 + 50.0**2 + 0.0**2) ** 2.0
+        kinEneNew = kinEneNew + particles["m"][k] * np.sqrt(11.0**2 + 11.0**2 + 0**2) ** 2 * 0.5
+        potEneNew = potEneNew + particles["m"][k] * 9.81 + 0.0
 
     # call function to be tested
-    particles = DFAfunC.updatePositionC(cfg['GENERAL'], particles, dem, force, fields, typeStop=typeStop)
-#    print('sph', particles['peakForceSPH'], sphForceNew)
-
-    assert np.array_equal(particles['m'], np.asarray([10., 10., 10.]))
-    assert np.array_equal(particles['x'], np.array([6., 7., 8.]))
-    assert np.array_equal(particles['y'], np.asarray([8., 9., 10.]))
-    assert np.array_equal(particles['z'], np.asarray([1., 1., 1.]))
-    assert np.allclose(particles['ux'], np.asarray([11., 11., 11.]), atol=1.e-4)
-    assert np.allclose(particles['uy'], np.asarray([11., 11., 11.]), atol=1.e-4)
-    assert np.allclose(particles['uAcc'], np.asarray([uAcc, uAcc, uAcc]), atol=1.e-4)
-    assert np.array_equal(particles['uz'], np.asarray([0., 0., 0.]))
-    assert particles['massEntrained'] == 0.0
-    assert particles['nPart'] == 3
-    assert (kinEneNew- 1.e-4) < particles['kineticEne'] < (kinEneNew+1.e-4)
-    assert (potEneNew-1.e-4) < particles['potentialEne'] < (potEneNew +1.e-4)
-    assert particles['iterate'] == True
+    particles = DFAfunC.updatePositionC(cfg["GENERAL"], particles, dem, force, fields, typeStop=typeStop)
+    #    print('sph', particles['peakForceSPH'], sphForceNew)
+
+    assert np.array_equal(particles["m"], np.asarray([10.0, 10.0, 10.0]))
+    assert np.array_equal(particles["x"], np.array([6.0, 7.0, 8.0]))
+    assert np.array_equal(particles["y"], np.asarray([8.0, 9.0, 10.0]))
+    assert np.array_equal(particles["z"], np.asarray([1.0, 1.0, 1.0]))
+    assert np.allclose(particles["ux"], np.asarray([11.0, 11.0, 11.0]), atol=1.0e-4)
+    assert np.allclose(particles["uy"], np.asarray([11.0, 11.0, 11.0]), atol=1.0e-4)
+    assert np.allclose(particles["uAcc"], np.asarray([uAcc, uAcc, uAcc]), atol=1.0e-4)
+    assert np.array_equal(particles["uz"], np.asarray([0.0, 0.0, 0.0]))
+    assert particles["massEntrained"] == 0.0
+    assert particles["nPart"] == 3
+    assert (kinEneNew - 1.0e-4) < particles["kineticEne"] < (kinEneNew + 1.0e-4)
+    assert (potEneNew - 1.0e-4) < particles["potentialEne"] < (potEneNew + 1.0e-4)
+    assert particles["iterate"] == True
 
 
 def test_computeTrajectoryAngle():
@@ -397,10 +490,10 @@ def test_computeTrajectoryAngle():
     zPartArray0 = np.array([10.0, 9.0, 8.0])
     s = np.array([10.0, 10.0, 0.0, 10.0])
     z = np.array([0.0, 0.0, 0.0, 1.0])
-    particles = {'nPart': nPart, 'parentID': parentID, 'trajectoryLengthXY': s, 'z': z}
+    particles = {"nPart": nPart, "parentID": parentID, "trajectoryLengthXY": s, "z": z}
     particles = DFAfunC.computeTrajectoryAngleC(particles, zPartArray0)
-#    print(particles['trajectoryAngle'])
-    gamma = particles['trajectoryAngle']
+    #    print(particles['trajectoryAngle'])
+    gamma = particles["trajectoryAngle"]
     assert gamma[2] == 0
     assert gamma[0] == 45
     assert gamma[1] == pytest.approx(41.9872125, rel=1e-6)
@@ -409,21 +502,21 @@ def test_computeTrajectoryAngle():
 
 def test_initializeBondsC():
     nPart = 3
-    x = np.array([0., 1., 0.])
-    y = np.array([0., 0., 1.])
-    z = np.array([0., 0., 0.])
+    x = np.array([0.0, 1.0, 0.0])
+    y = np.array([0.0, 0.0, 1.0])
+    z = np.array([0.0, 0.0, 0.0])
     # original triangulation
     triangles = tri.Triangulation(x, y)
-    particles = {'nPart': nPart, 'x': x, 'y': y, 'z': z}
+    particles = {"nPart": nPart, "x": x, "y": y, "z": z}
     particles = DFAfunC.initializeBondsC(particles, triangles)
-#    print(triangles.triangles)
-#    print(triangles.edges)
-#    print(particles['bondStart'])
-#    print(particles['bondDist'])
-#    print(particles['bondPart'])
-    bondStart = particles['bondStart']
-    bondDist = particles['bondDist']
-    bondPart = particles['bondPart']
+    #    print(triangles.triangles)
+    #    print(triangles.edges)
+    #    print(particles['bondStart'])
+    #    print(particles['bondDist'])
+    #    print(particles['bondPart'])
+    bondStart = particles["bondStart"]
+    bondDist = particles["bondDist"]
+    bondPart = particles["bondPart"]
     assert np.array_equal(bondStart, np.asarray([0, 2, 4, 6]))
     for k in range(nPart):
         # loop on all bonded particles
@@ -445,30 +538,30 @@ def test_initializeBondsC():
 
 def test_removeBondsC():
     nPart = 3
-    x = np.array([0., 1., 0.])
-    y = np.array([0., 0., 1.])
-    z = np.array([0., 0., 0.])
+    x = np.array([0.0, 1.0, 0.0])
+    y = np.array([0.0, 0.0, 1.0])
+    z = np.array([0.0, 0.0, 0.0])
     # original triangulation
     triangles = tri.Triangulation(x, y)
-    particles = {'nPart': nPart, 'x': x, 'y': y, 'z': z}
+    particles = {"nPart": nPart, "x": x, "y": y, "z": z}
     particles = DFAfunC.initializeBondsC(particles, triangles)
-#    print(particles['bondStart'])
-#    print(particles['bondDist'])
-#    print(particles['bondPart'])
+    #    print(particles['bondStart'])
+    #    print(particles['bondDist'])
+    #    print(particles['bondPart'])
 
     # now remove one particle (here particle 1):
-    keepParticle = np.array([1., 0., 1.])
+    keepParticle = np.array([1.0, 0.0, 1.0])
     nRemove = 1
     nBondRemove = DFAfunC.countRemovedBonds(particles, keepParticle, nRemove)
-#    print(nBondRemove)
+    #    print(nBondRemove)
     assert nBondRemove == 4
     particles = DFAfunC.removedBonds(particles, keepParticle, nRemove, nBondRemove)
-#    print(particles['bondStart'])
-#    print(particles['bondDist'])
-#    print(particles['bondPart'])
-    bondStart = particles['bondStart']
-    bondDist = particles['bondDist']
-    bondPart = particles['bondPart']
+    #    print(particles['bondStart'])
+    #    print(particles['bondDist'])
+    #    print(particles['bondPart'])
+    bondStart = particles["bondStart"]
+    bondDist = particles["bondDist"]
+    bondPart = particles["bondPart"]
     assert np.array_equal(bondStart, np.asarray([0, 1, 1, 2]))
     for k in range(nPart):
         # loop on all bonded particles
@@ -490,38 +583,53 @@ def test_removeBondsC():
 
 def test_computeCohesionForceC():
     cfg = configparser.ConfigParser()
-    cfg['GENERAL'] = {'cohesiveSurfaceTension': '50000', 'cohesionMaxStrain': '0.2', 'minDistCohesion': '1.0e-3'}
+    cfg["GENERAL"] = {
+        "cohesiveSurfaceTension": "50000",
+        "cohesionMaxStrain": "0.2",
+        "minDistCohesion": "1.0e-3",
+    }
     nPart = 3
-    x = np.array([0., 1., 0.])
-    y = np.array([0., 0., 1.])
-    z = np.array([0., 0., 0.])
-    ux = np.array([0., 0., 0.])
-    uy = np.array([0., 0., 0.])
-    uz = np.array([0., 0., 0.])
-    m = np.array([1., 1., 1.])
-    h = np.array([1., 1., 1.])
+    x = np.array([0.0, 1.0, 0.0])
+    y = np.array([0.0, 0.0, 1.0])
+    z = np.array([0.0, 0.0, 0.0])
+    ux = np.array([0.0, 0.0, 0.0])
+    uy = np.array([0.0, 0.0, 0.0])
+    uz = np.array([0.0, 0.0, 0.0])
+    m = np.array([1.0, 1.0, 1.0])
+    h = np.array([1.0, 1.0, 1.0])
     force = {}
-    force['forceSPHX'] = np.zeros(np.shape(x))
-    force['forceSPHY'] = np.zeros(np.shape(x))
-    force['forceSPHZ'] = np.zeros(np.shape(x))
+    force["forceSPHX"] = np.zeros(np.shape(x))
+    force["forceSPHY"] = np.zeros(np.shape(x))
+    force["forceSPHZ"] = np.zeros(np.shape(x))
     # original triangulation
     triangles = tri.Triangulation(x, y)
-    particles = {'nPart': nPart, 'x': x, 'y': y, 'z': z, 'ux': ux, 'uy': uy, 'uz': uz, 'm': m, 'h': h, 'dt': 0.05}
+    particles = {
+        "nPart": nPart,
+        "x": x,
+        "y": y,
+        "z": z,
+        "ux": ux,
+        "uy": uy,
+        "uz": uz,
+        "m": m,
+        "h": h,
+        "dt": 0.05,
+    }
     particles = DFAfunC.initializeBondsC(particles, triangles)
     # bond breaking
-    particles['x'][1] = 1.21
-#    print('Here')
-    force, particles = DFAfunC.computeCohesionForceC(cfg['GENERAL'], particles, force)
-#    print('Here')
-#    print(particles['bondStart'])
-#    print(particles['bondDist'])
-#    print(particles['bondPart'])
-#    print(force['forceSPHX'])
-#    print(force['forceSPHY'])
-#    print(force['forceSPHZ'])
-    bondStart = particles['bondStart']
-    bondDist = particles['bondDist']
-    bondPart = particles['bondPart']
+    particles["x"][1] = 1.21
+    #    print('Here')
+    force, particles = DFAfunC.computeCohesionForceC(cfg["GENERAL"], particles, force)
+    #    print('Here')
+    #    print(particles['bondStart'])
+    #    print(particles['bondDist'])
+    #    print(particles['bondPart'])
+    #    print(force['forceSPHX'])
+    #    print(force['forceSPHY'])
+    #    print(force['forceSPHZ'])
+    bondStart = particles["bondStart"]
+    bondDist = particles["bondDist"]
+    bondPart = particles["bondPart"]
     for k in range(nPart):
         # loop on all bonded particles
         neighbors = list()
@@ -544,30 +652,41 @@ def test_computeCohesionForceC():
             assert neighbors == [0, 1]
     bondDist.sort()
     assert np.array_equal(bondDist, np.asarray([-1, -1, 1, 1, math.sqrt(2), math.sqrt(2)]))
-    assert force['forceSPHX'][0] == 0
-    assert force['forceSPHY'][0] == 0
-    assert force['forceSPHX'][1] < 0
-    assert force['forceSPHY'][1] > 0
-    assert force['forceSPHX'][2] > 0
-    assert force['forceSPHY'][2] < 0
-    assert force['forceSPHX'][1] == -force['forceSPHX'][2]
-    assert force['forceSPHY'][1] == -force['forceSPHY'][2]
+    assert force["forceSPHX"][0] == 0
+    assert force["forceSPHY"][0] == 0
+    assert force["forceSPHX"][1] < 0
+    assert force["forceSPHY"][1] > 0
+    assert force["forceSPHX"][2] > 0
+    assert force["forceSPHY"][2] < 0
+    assert force["forceSPHX"][1] == -force["forceSPHX"][2]
+    assert force["forceSPHY"][1] == -force["forceSPHY"][2]
 
     # bound not breaking
-    particles['x'][1] = 1
-    particles = {'nPart': nPart, 'x': x, 'y': y, 'z': z, 'ux': ux, 'uy': uy, 'uz': uz, 'm': m, 'h': h, 'dt': 0.05}
+    particles["x"][1] = 1
+    particles = {
+        "nPart": nPart,
+        "x": x,
+        "y": y,
+        "z": z,
+        "ux": ux,
+        "uy": uy,
+        "uz": uz,
+        "m": m,
+        "h": h,
+        "dt": 0.05,
+    }
     particles = DFAfunC.initializeBondsC(particles, triangles)
-    particles['x'][1] = 1.09
-    force, particles = DFAfunC.computeCohesionForceC(cfg['GENERAL'], particles, force)
-#    print(particles['bondStart'])
-#    print(particles['bondDist'])
-#    print(particles['bondPart'])
-#    print(force['forceSPHX'])
-#    print(force['forceSPHY'])
-#    print(force['forceSPHZ'])
-    bondStart = particles['bondStart']
-    bondDist = particles['bondDist']
-    bondPart = particles['bondPart']
+    particles["x"][1] = 1.09
+    force, particles = DFAfunC.computeCohesionForceC(cfg["GENERAL"], particles, force)
+    #    print(particles['bondStart'])
+    #    print(particles['bondDist'])
+    #    print(particles['bondPart'])
+    #    print(force['forceSPHX'])
+    #    print(force['forceSPHY'])
+    #    print(force['forceSPHZ'])
+    bondStart = particles["bondStart"]
+    bondDist = particles["bondDist"]
+    bondPart = particles["bondPart"]
     for k in range(nPart):
         # loop on all bonded particles
         neighbors = list()
@@ -585,15 +704,16 @@ def test_computeCohesionForceC():
     bondDist.sort()
     assert np.array_equal(bondDist, np.asarray([1, 1, 1, 1, math.sqrt(2), math.sqrt(2)]))
 
-    assert force['forceSPHX'][0] > 0
-    assert force['forceSPHY'][0] == 0
-    assert force['forceSPHX'][1] < 0
-    assert force['forceSPHY'][1] > 0
-    assert force['forceSPHX'][2] > 0
-    assert force['forceSPHY'][2] < 0
-    assert force['forceSPHY'][1] == -force['forceSPHY'][2]
+    assert force["forceSPHX"][0] > 0
+    assert force["forceSPHY"][0] == 0
+    assert force["forceSPHX"][1] < 0
+    assert force["forceSPHY"][1] > 0
+    assert force["forceSPHX"][2] > 0
+    assert force["forceSPHY"][2] < 0
+    assert force["forceSPHY"][1] == -force["forceSPHY"][2]
+
 
-'''
+"""
 TODO: When calling pytest, the following function raises an error ("Fatal Python error: Aborted")
 (see issue #1002?)
 
@@ -652,4 +772,4 @@ def test_updateFieldsC():
     atol = 1e-10
     assert np.allclose(fields['dmDet'], dmDet_calculated2, atol=atol)
 #    print(fields['dmDet'])
-    '''
+    """
diff --git a/avaframe/tests/test_com1DFA.py b/avaframe/tests/test_com1DFA.py
index e485ee37f..36090bfe8 100644
--- a/avaframe/tests/test_com1DFA.py
+++ b/avaframe/tests/test_com1DFA.py
@@ -1,5 +1,5 @@
 """
-    Pytest for module com1DFA
+Pytest for module com1DFA
 """
 
 import configparser
@@ -38,7 +38,11 @@ def test_prepareInputData(tmp_path):
     inputSimFiles["muFile"] = None
     inputSimFiles["xiFile"] = None
     cfg = configparser.ConfigParser()
-    cfg["GENERAL"] = {"secRelArea": "False", "simTypeActual": "ent", "avalancheDir": str(avaDir)}
+    cfg["GENERAL"] = {
+        "secRelArea": "False",
+        "simTypeActual": "ent",
+        "avalancheDir": str(avaDir),
+    }
     cfg["GENERAL"]["relThFromFile"] = "False"
     cfg["INPUT"] = {"DEM": "avaAlr.tif"}
     cfg["INPUT"]["relThFile"] = ""
@@ -285,7 +289,11 @@ def test_prepareInputData(tmp_path):
     inputSimFiles["muFile"] = None
     inputSimFiles["xiFile"] = None
     cfg = configparser.ConfigParser()
-    cfg["GENERAL"] = {"secRelArea": "False", "simTypeActual": "null", "avalancheDir": str(avaDir)}
+    cfg["GENERAL"] = {
+        "secRelArea": "False",
+        "simTypeActual": "null",
+        "avalancheDir": str(avaDir),
+    }
     cfg["INPUT"] = {"DEM": "avaAlr.tif"}
     cfg["INPUT"]["relThFile"] = ""
 
@@ -309,11 +317,19 @@ def test_prepareReleaseEntrainment(tmp_path):
         "secondaryRelThPercentVariation": "0.7",
         "simTypeActual": "null",
     }
-    cfg["INPUT"] = {"secondaryRelThThickness": "1.2523", "secondaryRelThId": "0", "thFromIni": ""}
+    cfg["INPUT"] = {
+        "secondaryRelThThickness": "1.2523",
+        "secondaryRelThId": "0",
+        "thFromIni": "",
+    }
 
     inputSimLines = {}
     inputSimLines["entResInfo"] = {"flagSecondaryRelease": "Yes", "flagEnt": "No"}
-    inputSimLines["releaseLine"] = {"thickness": ["None", "None"], "type": "Release", "id": ["0", "1"]}
+    inputSimLines["releaseLine"] = {
+        "thickness": ["None", "None"],
+        "type": "Release",
+        "id": ["0", "1"],
+    }
     inputSimLines["relThField"] = ""
     inputSimLines["secondaryReleaseLine"] = {
         "thickness": ["1.2523"],
@@ -346,8 +362,16 @@ def test_prepareReleaseEntrainment(tmp_path):
 
     inputSimLines = {}
     inputSimLines["entResInfo"] = {"flagSecondaryRelease": "Yes", "flagEnt": "No"}
-    inputSimLines["releaseLine"] = {"thickness": ["1.78", "4.328"], "type": "release", "id": ["0", "1"]}
-    inputSimLines["secondaryReleaseLine"] = {"thickness": ["None"], "type": "Secondary release", "id": ["0"]}
+    inputSimLines["releaseLine"] = {
+        "thickness": ["1.78", "4.328"],
+        "type": "release",
+        "id": ["0", "1"],
+    }
+    inputSimLines["secondaryReleaseLine"] = {
+        "thickness": ["None"],
+        "type": "Secondary release",
+        "id": ["0"],
+    }
     inputSimLines["relThField"] = ""
     rel = pathlib.Path(tmp_path, "release1PF_test.shp")
 
@@ -379,8 +403,16 @@ def test_prepareReleaseEntrainment(tmp_path):
 
     inputSimLines = {}
     inputSimLines["entResInfo"] = {"flagSecondaryRelease": "Yes", "flagEnt": "No"}
-    inputSimLines["releaseLine"] = {"thickness": ["1.", "2."], "type": "release", "id": ["0", "1"]}
-    inputSimLines["secondaryReleaseLine"] = {"thickness": ["2.7"], "type": "Secondary release", "id": ["0"]}
+    inputSimLines["releaseLine"] = {
+        "thickness": ["1.", "2."],
+        "type": "release",
+        "id": ["0", "1"],
+    }
+    inputSimLines["secondaryReleaseLine"] = {
+        "thickness": ["2.7"],
+        "type": "Secondary release",
+        "id": ["0"],
+    }
     rel = pathlib.Path(tmp_path, "release1PF_test.shp")
 
     # call function to be tested
@@ -415,7 +447,11 @@ def test_prepareReleaseEntrainment(tmp_path):
     # setup required inputs
     inputSimLines = {}
     inputSimLines["entResInfo"] = {"flagSecondaryRelease": "No", "flagEnt": "No"}
-    inputSimLines["releaseLine"] = {"thickness": ["1.78", "4.328"], "type": "release", "id": ["0", "1"]}
+    inputSimLines["releaseLine"] = {
+        "thickness": ["1.78", "4.328"],
+        "type": "release",
+        "id": ["0", "1"],
+    }
     rel = pathlib.Path(tmp_path, "release1PF_test.shp")
     cfg["GENERAL"]["relThFromShp"] = "False"
     cfg["GENERAL"]["relTh"] = "1.32"
@@ -445,8 +481,16 @@ def test_prepareReleaseEntrainment(tmp_path):
     }
     inputSimLines = {}
     inputSimLines["entResInfo"] = {"flagSecondaryRelease": "No", "flagEnt": "Yes"}
-    inputSimLines["releaseLine"] = {"thickness": ["None", "None"], "type": "Release", "id": ["0", "1"]}
-    inputSimLines["entLine"] = {"thickness": ["1.20", "0.9"], "type": "Entrainment", "id": ["0", "1"]}
+    inputSimLines["releaseLine"] = {
+        "thickness": ["None", "None"],
+        "type": "Release",
+        "id": ["0", "1"],
+    }
+    inputSimLines["entLine"] = {
+        "thickness": ["1.20", "0.9"],
+        "type": "Entrainment",
+        "id": ["0", "1"],
+    }
     relName5, inputSimLines5, badName5 = com1DFA.prepareReleaseEntrainment(cfg, rel, inputSimLines)
 
     assert relName5 == "release1PF_test"
@@ -463,7 +507,11 @@ def test_setThickness():
 
     # setup required input
     cfg = configparser.ConfigParser()
-    cfg["GENERAL"] = {"entThFromShp": "False", "entTh": "1.0", "entThPercentVariation": ""}
+    cfg["GENERAL"] = {
+        "entThFromShp": "False",
+        "entTh": "1.0",
+        "entThPercentVariation": "",
+    }
     cfg["INPUT"] = {"thFromIni": ""}
 
     lineTh = {
@@ -646,7 +694,12 @@ def test_initializeMassEnt():
     dem["header"]["yllcenter"] = 0.0
 
     cfg = configparser.ConfigParser()
-    cfg["GENERAL"] = {"rhoEnt": "200.0", "entTempRef": "-10.", "cpIce": "2050.", "TIni": "-10."}
+    cfg["GENERAL"] = {
+        "rhoEnt": "200.0",
+        "entTempRef": "-10.",
+        "cpIce": "2050.",
+        "TIni": "-10.",
+    }
 
     simTypeActual = "entres"
     dirName = pathlib.Path(__file__).parents[0]
@@ -669,7 +722,12 @@ def test_initializeMassEnt():
 
     # call function to be tested
     entrMassRaster, entrEnthRaster, reportAreaInfo = com1DFA.initializeMassEnt(
-        dem, simTypeActual, entLine, reportAreaInfo, thresholdPointInPoly, cfg["GENERAL"]
+        dem,
+        simTypeActual,
+        entLine,
+        reportAreaInfo,
+        thresholdPointInPoly,
+        cfg["GENERAL"],
     )
     testData = np.zeros((nrows, ncols))
     testData[0:11, 0:11] = 1.0 * 200.0
@@ -687,7 +745,12 @@ def test_initializeMassEnt():
     # call function to be tested
     simTypeActual = "res"
     entrMassRaster, entrEnthRaster, reportAreaInfo = com1DFA.initializeMassEnt(
-        dem, simTypeActual, entLine, reportAreaInfo, thresholdPointInPoly, cfg["GENERAL"]
+        dem,
+        simTypeActual,
+        entLine,
+        reportAreaInfo,
+        thresholdPointInPoly,
+        cfg["GENERAL"],
     )
 
     assert np.array_equal(entrMassRaster, np.zeros((nrows, ncols)))
@@ -700,11 +763,14 @@ def test_initializeResistance():
     # setup required input
     cfg = configparser.ConfigParser()
     cfg["GENERAL"] = {
-        "cRes": 0.003,
-        "ResistanceModel": "cRes",
+        "cResH": 0.003,
+        "ResistanceModel": "default",
         "detK": 10,
         "detrainment": False,
-        "detWithoutRes": False,
+        "forestVMin": 6.0,
+        "forestThMin": 0.6,
+        "forestVMax": 40.0,
+        "forestThMax": 10.0,
     }
 
     nrows = 11
@@ -737,7 +803,12 @@ def test_initializeResistance():
     dem["header"]["xllcenter"] = 0.0
     dem["header"]["yllcenter"] = 0.0
     cResRaster, detRaster, reportAreaInfo = com1DFA.initializeResistance(
-        cfg["GENERAL"], dem, simTypeActual, resLine, reportAreaInfo, thresholdPointInPoly
+        cfg["GENERAL"],
+        dem,
+        simTypeActual,
+        resLine,
+        reportAreaInfo,
+        thresholdPointInPoly,
     )
     testArray = np.zeros((nrows, ncols))
     testArray[0:11, 0:11] = 0.003
@@ -753,14 +824,22 @@ def test_initializeResistance():
     assert reportAreaInfo["detrainment"] == "No"
 
     cfg["GENERAL"] = {
-        "cRes": 0.003,
-        "ResistanceModel": "cRes",
+        "cResH": 0.003,
+        "ResistanceModel": "default",
         "detK": 10.0,
         "detrainment": True,
-        "detWithoutRes": False,
+        "forestVMin": 6.0,
+        "forestThMin": 0.6,
+        "forestVMax": 40.0,
+        "forestThMax": 10.0,
     }
     cResRaster, detRaster, reportAreaInfo = com1DFA.initializeResistance(
-        cfg["GENERAL"], dem, simTypeActual, resLine, reportAreaInfo, thresholdPointInPoly
+        cfg["GENERAL"],
+        dem,
+        simTypeActual,
+        resLine,
+        reportAreaInfo,
+        thresholdPointInPoly,
     )
     detTestArray = np.zeros((nrows, ncols))
     detTestArray[0:11, 0:11] = 10.0
@@ -773,38 +852,26 @@ def test_initializeResistance():
     assert reportAreaInfo["detrainment"] == "Yes"
 
     cfg["GENERAL"] = {
-        "cRes": 0.003,
+        "cResH": 0.003,
         "ResistanceModel": "cRes",
         "detK": 10.0,
         "detrainment": True,
-        "detWithoutRes": True,
+        "forestVMin": 6.0,
+        "forestThMin": 0.6,
+        "forestVMax": 40.0,
+        "forestThMax": 10.0,
     }
-    cResRaster, detRaster, reportAreaInfo = com1DFA.initializeResistance(
-        cfg["GENERAL"], dem, simTypeActual, resLine, reportAreaInfo, thresholdPointInPoly
-    )
-    detTestArray = np.zeros((nrows, ncols))
-    detTestArray[0:11, 0:11] = 10.0
 
-    assert np.array_equal(cResRaster, np.zeros((nrows, ncols)))
-    assert np.array_equal(detRaster, detTestArray)
-    assert np.sum(detRaster) == 1210.0
-    assert np.sum(cResRaster) == 0.0
-    assert reportAreaInfo["resistance"] == "No"
-    assert reportAreaInfo["detrainment"] == "Yes"
-
-    cfg["GENERAL"] = {
-        "cResH": 0.003,
-        "ResistanceModel": "cResH",
-        "detK": 10.0,
-        "detrainment": False,
-        "detWithoutRes": False,
-    }
-    cResRaster, detRaster, reportAreaInfo = com1DFA.initializeResistance(
-        cfg["GENERAL"], dem, simTypeActual, resLine, reportAreaInfo, thresholdPointInPoly
-    )
-    assert np.array_equal(cResRaster, testArray)
-    assert np.sum(cResRaster) == 0.363
-    assert reportAreaInfo["resistance"] == "Yes"
+    with pytest.raises(AssertionError) as e:
+        assert com1DFA.initializeResistance(
+            cfg["GENERAL"],
+            dem,
+            simTypeActual,
+            resLine,
+            reportAreaInfo,
+            thresholdPointInPoly,
+        )
+    assert "Resistance model cres not a valid option" in str(e.value)
 
 
 def test_setDEMOriginToZero():
@@ -848,7 +915,11 @@ def test_initializeMesh():
 
     demOri = {"header": demHeader, "rasterData": demData}
     cfg = configparser.ConfigParser()
-    cfg["GENERAL"] = {"sphKernelRadius": "0.5", "meshCellSizeThreshold": "0.0001", "meshCellSize": "1."}
+    cfg["GENERAL"] = {
+        "sphKernelRadius": "0.5",
+        "meshCellSizeThreshold": "0.0001",
+        "meshCellSize": "1.",
+    }
     num = 1
 
     # setup testResults
@@ -1015,7 +1086,12 @@ def test_appendFieldsParticles():
         "y": np.asarray([0.0, 5.0, 0.0]),
         "m": np.asarray([0.0, 4.0, 0.0]),
     }
-    fields = {"ppr": np.ones((3, 3)), "pft": np.ones((3, 3)), "pfv": np.ones((3, 3)), "FT": np.ones((3, 3))}
+    fields = {
+        "ppr": np.ones((3, 3)),
+        "pft": np.ones((3, 3)),
+        "pfv": np.ones((3, 3)),
+        "FT": np.ones((3, 3)),
+    }
     resTypes = ["ppr", "pfv", "pft", "particles"]
 
     # call function to be tested
@@ -1076,9 +1152,63 @@ def test_releaseSecRelArea():
     secRelRaster1 = np.zeros((demHeader["nrows"], demHeader["ncols"]))
     secRelRaster1[1, 1] = 0.5
     secondaryReleaseInfo = {
-        "x": np.asarray([1.5, 2.5, 2.5, 1.5, 1.5, 7.4, 8.5, 8.5, 7.4, 7.4, 9.5, 10.5, 10.5, 9.5, 9.5]),
-        "y": np.asarray([1.5, 1.5, 2.5, 2.5, 1.5, 7.4, 7.4, 8.5, 8.5, 7.4, 9.5, 9.5, 10.5, 10.5, 9.5]),
-        "z": np.asarray([1.5, 1.5, 2.5, 2.5, 1.5, 7.4, 7.4, 8.5, 8.5, 7.4, 9.5, 9.5, 10.5, 10.5, 9.5]),
+        "x": np.asarray(
+            [
+                1.5,
+                2.5,
+                2.5,
+                1.5,
+                1.5,
+                7.4,
+                8.5,
+                8.5,
+                7.4,
+                7.4,
+                9.5,
+                10.5,
+                10.5,
+                9.5,
+                9.5,
+            ]
+        ),
+        "y": np.asarray(
+            [
+                1.5,
+                1.5,
+                2.5,
+                2.5,
+                1.5,
+                7.4,
+                7.4,
+                8.5,
+                8.5,
+                7.4,
+                9.5,
+                9.5,
+                10.5,
+                10.5,
+                9.5,
+            ]
+        ),
+        "z": np.asarray(
+            [
+                1.5,
+                1.5,
+                2.5,
+                2.5,
+                1.5,
+                7.4,
+                7.4,
+                8.5,
+                8.5,
+                7.4,
+                9.5,
+                9.5,
+                10.5,
+                10.5,
+                9.5,
+            ]
+        ),
         "Start": np.asarray([0, 5, 10]),
         "Length": np.asarray([5, 5, 5]),
         "Name": ["secRel1", "secRel2", "secRel3"],
@@ -1140,14 +1270,17 @@ def test_releaseSecRelArea():
     assert particles["mTot"] == 2700.0
     assert particles2["nPart"] == 11
     assert np.array_equal(
-        particles2["x"], np.asarray([6.0, 7.0, 9.1, 6.75, 7.25, 6.75, 7.25, 8.75, 9.25, 8.75, 9.25])
+        particles2["x"],
+        np.asarray([6.0, 7.0, 9.1, 6.75, 7.25, 6.75, 7.25, 8.75, 9.25, 8.75, 9.25]),
     )
     assert np.array_equal(
-        particles2["y"], np.asarray([6.0, 7.0, 9.1, 6.75, 6.75, 7.25, 7.25, 8.75, 8.75, 9.25, 9.25])
+        particles2["y"],
+        np.asarray([6.0, 7.0, 9.1, 6.75, 6.75, 7.25, 7.25, 8.75, 8.75, 9.25, 9.25]),
     )
     assert np.array_equal(zPartArray0New2, np.asarray([1, 2, 3, 1, 1, 1, 1, 1, 1, 1, 1]))
     assert np.array_equal(
-        particles2["m"], np.asarray([1250.0, 1250.0, 1250.0, 50.0, 50.0, 50.0, 50.0, 25.0, 25.0, 25.0, 25.0])
+        particles2["m"],
+        np.asarray([1250.0, 1250.0, 1250.0, 50.0, 50.0, 50.0, 50.0, 25.0, 25.0, 25.0, 25.0]),
     )
     assert particles2["mTot"] == 4050.0
 
@@ -1286,6 +1419,7 @@ def test_initializeParticles():
         "totalEnthalpy",
         "velocityMag",
         "nExitedParticles",
+        "tPlot",
     ]
 
     # call function to be tested
@@ -1308,9 +1442,18 @@ def test_initializeParticles():
     assert all(key in particles for key in dictKeys)
 
     assert particles["nPart"] == 9
-    assert np.array_equal(particles["x"], np.asarray([6.25, 6.75, 7.25, 6.25, 6.25, 6.75, 7.25, 6.75, 7.25]))
-    assert np.array_equal(particles["y"], np.asarray([6.25, 6.25, 6.25, 6.75, 7.25, 6.75, 6.75, 7.25, 7.25]))
-    assert np.array_equal(particles["m"], np.asarray([50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0]))
+    assert np.array_equal(
+        particles["x"],
+        np.asarray([6.25, 6.75, 7.25, 6.25, 6.25, 6.75, 7.25, 6.75, 7.25]),
+    )
+    assert np.array_equal(
+        particles["y"],
+        np.asarray([6.25, 6.25, 6.25, 6.75, 7.25, 6.75, 6.75, 7.25, 7.25]),
+    )
+    assert np.array_equal(
+        particles["m"],
+        np.asarray([50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0]),
+    )
     assert particles["mTot"] == 450.0
     assert np.sum(particles["ux"]) == 0.0
 
@@ -1336,9 +1479,18 @@ def test_initializeParticles():
     #    print("particles", particles)
 
     assert particles["nPart"] == 9
-    assert np.array_equal(particles["x"], np.asarray([6.25, 6.75, 7.25, 6.25, 6.25, 6.75, 7.25, 6.75, 7.25]))
-    assert np.array_equal(particles["y"], np.asarray([6.25, 6.25, 6.25, 6.75, 7.25, 6.75, 6.75, 7.25, 7.25]))
-    assert np.array_equal(particles["m"], np.asarray([50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0]))
+    assert np.array_equal(
+        particles["x"],
+        np.asarray([6.25, 6.75, 7.25, 6.25, 6.25, 6.75, 7.25, 6.75, 7.25]),
+    )
+    assert np.array_equal(
+        particles["y"],
+        np.asarray([6.25, 6.25, 6.25, 6.75, 7.25, 6.75, 6.75, 7.25, 7.25]),
+    )
+    assert np.array_equal(
+        particles["m"],
+        np.asarray([50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0]),
+    )
     assert particles["mTot"] == 450.0
     assert np.sum(particles["ux"]) == 0.0
 
@@ -1367,9 +1519,18 @@ def test_initializeParticles():
     #    print("particles", particles)
 
     assert particles["nPart"] == 9
-    assert np.array_equal(particles["x"], np.asarray([6.25, 6.75, 7.25, 6.25, 6.25, 6.75, 7.25, 6.75, 7.25]))
-    assert np.array_equal(particles["y"], np.asarray([6.25, 6.25, 6.25, 6.75, 7.25, 6.75, 6.75, 7.25, 7.25]))
-    assert np.array_equal(particles["m"], np.asarray([50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0]))
+    assert np.array_equal(
+        particles["x"],
+        np.asarray([6.25, 6.75, 7.25, 6.25, 6.25, 6.75, 7.25, 6.75, 7.25]),
+    )
+    assert np.array_equal(
+        particles["y"],
+        np.asarray([6.25, 6.25, 6.25, 6.75, 7.25, 6.75, 6.75, 7.25, 7.25]),
+    )
+    assert np.array_equal(
+        particles["m"],
+        np.asarray([50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0]),
+    )
     assert particles["mTot"] == 450.0
     assert np.sum(particles["ux"]) == 0.0
 
@@ -1382,6 +1543,7 @@ def test_writeMBFile(tmp_path):
     infoDict["massEntrained"] = np.asarray([0, 0, 10, 20, 30])
     infoDict["massDetrained"] = np.asarray([0, 0, 0, 0, 0])
     infoDict["massTotal"] = np.asarray([60.0, 60.0, 70.0, 90.0, 120.0])
+    infoDict["pfvTimeMax"] = np.asarray([0, 0, 0, 0, 0])
     avaName = "data/avaTest"
     avaDir = pathlib.Path(tmp_path, avaName)
     logName = "simTestName"
@@ -1399,7 +1561,7 @@ def test_writeMBFile(tmp_path):
     assert np.array_equal(mbInfo[:, 3], infoDict["massDetrained"])
     assert np.array_equal(mbInfo[:, 1], infoDict["massTotal"])
     assert mbInfo.shape[0] == 5
-    assert mbInfo.shape[1] == 4
+    assert mbInfo.shape[1] == 5
 
     infoDict["massEntrained"] = np.asarray([0, 0, 0, 0, 0])
     infoDict["massDetrained"] = np.asarray([0, 10, 0, 30, 0])
@@ -1414,7 +1576,7 @@ def test_writeMBFile(tmp_path):
     assert np.array_equal(mbInfo[:, 3], infoDict["massDetrained"])
     assert np.array_equal(mbInfo[:, 1], infoDict["massTotal"])
     assert mbInfo.shape[0] == 5
-    assert mbInfo.shape[1] == 4
+    assert mbInfo.shape[1] == 5
 
     infoDict["massEntrained"] = np.asarray([0, 20, 0, 0, 10])
     infoDict["massDetrained"] = np.asarray([0, 10, 0, 30, 0])
@@ -1429,7 +1591,7 @@ def test_writeMBFile(tmp_path):
     assert np.array_equal(mbInfo[:, 3], infoDict["massDetrained"])
     assert np.array_equal(mbInfo[:, 1], infoDict["massTotal"])
     assert mbInfo.shape[0] == 5
-    assert mbInfo.shape[1] == 4
+    assert mbInfo.shape[1] == 5
 
 
 def test_savePartToPickle(tmp_path):
@@ -1523,17 +1685,46 @@ def test_exportFields(tmp_path):
     pft = np.zeros((5, 5))
     pfv = np.zeros((5, 5))
     ppr = np.zeros((5, 5))
-    fields1 = {"ppr": ppr + 1, "pft": pft + 1, "pfv": pfv + 1, "FT": FT + 1, "pke": pke + 1}
-    fields2 = {"ppr": ppr + 2, "pft": pft + 2, "pfv": pfv + 2, "FT": FT + 2, "pke": pke + 2}
-    fields3 = {"ppr": ppr + 4, "pft": pft + 4, "pfv": pfv + 4, "FT": FT + 4, "pke": pke + 4}
-    fields4 = {"ppr": ppr + 5, "pft": pft + 5, "pfv": pfv + 5, "FT": FT + 5, "pke": pke + 5}
-    fields5 = {"ppr": ppr + 6, "pft": pft + 6, "pfv": pfv + 6, "FT": FT + 6, "pke": pke + 6}
+    fields1 = {
+        "ppr": ppr + 1,
+        "pft": pft + 1,
+        "pfv": pfv + 1,
+        "FT": FT + 1,
+        "pke": pke + 1,
+    }
+    fields2 = {
+        "ppr": ppr + 2,
+        "pft": pft + 2,
+        "pfv": pfv + 2,
+        "FT": FT + 2,
+        "pke": pke + 2,
+    }
+    fields3 = {
+        "ppr": ppr + 4,
+        "pft": pft + 4,
+        "pfv": pfv + 4,
+        "FT": FT + 4,
+        "pke": pke + 4,
+    }
+    fields4 = {
+        "ppr": ppr + 5,
+        "pft": pft + 5,
+        "pfv": pfv + 5,
+        "FT": FT + 5,
+        "pke": pke + 5,
+    }
+    fields5 = {
+        "ppr": ppr + 6,
+        "pft": pft + 6,
+        "pfv": pfv + 6,
+        "FT": FT + 6,
+        "pke": pke + 6,
+    }
     fieldsList = [fields1, fields2, fields3, fields4, fields5]
 
     # call function to be tested
-    com1DFA.exportFields(cfg, 10.00, fields2, dem, outDir, logName, TSave='intermediate')
-    com1DFA.exportFields(cfg, 40.00, fields5, dem, outDir, logName, TSave='final')
-
+    com1DFA.exportFields(cfg, 10.00, fields2, dem, outDir, logName, TSave="intermediate")
+    com1DFA.exportFields(cfg, 40.00, fields5, dem, outDir, logName, TSave="final")
 
     # read fields
     fieldDir = outDir / "peakFiles"
@@ -1565,11 +1756,11 @@ def test_exportFields(tmp_path):
     cfg["GENERAL"]["resType"] = ""
     cfg["REPORT"] = {"plotFields": "ppr|pft|pfv"}
 
-    com1DFA.exportFields(cfg, 0.00, fields1, dem, outDir2, logName, TSave='initial')
-    com1DFA.exportFields(cfg, 10.00, fields2, dem, outDir2, logName, TSave='intermediate')
-    com1DFA.exportFields(cfg, 15.00, fields3, dem, outDir2, logName, TSave='intermediate')
-    com1DFA.exportFields(cfg, 25.00, fields4, dem, outDir2, logName, TSave='intermediate')
-    com1DFA.exportFields(cfg, 40.00, fields5, dem, outDir2, logName, TSave='final')
+    com1DFA.exportFields(cfg, 0.00, fields1, dem, outDir2, logName, TSave="initial")
+    com1DFA.exportFields(cfg, 10.00, fields2, dem, outDir2, logName, TSave="intermediate")
+    com1DFA.exportFields(cfg, 15.00, fields3, dem, outDir2, logName, TSave="intermediate")
+    com1DFA.exportFields(cfg, 25.00, fields4, dem, outDir2, logName, TSave="intermediate")
+    com1DFA.exportFields(cfg, 40.00, fields5, dem, outDir2, logName, TSave="final")
 
     # read fields
     fieldDir = outDir2 / "peakFiles"
@@ -1595,7 +1786,11 @@ def test_initializeFields():
     # setup required inputs
     demHeader = {"nrows": 11, "ncols": 12, "cellsize": 1}
     areaRaster = np.ones((11, 12))
-    dem = {"header": demHeader, "headerNeighbourGrid": demHeader, "areaRaster": areaRaster}
+    dem = {
+        "header": demHeader,
+        "headerNeighbourGrid": demHeader,
+        "areaRaster": areaRaster,
+    }
     particles = {
         "x": np.asarray([1.0, 2.0, 3.0]),
         "y": np.asarray([1.0, 2.0, 3.0]),
@@ -1707,11 +1902,14 @@ def test_prepareVarSimDict(tmp_path, caplog):
         "entResInfo": {"flagEnt": "Yes", "flagRes": "Yes"},
         "demFile": avaDEM,
         "damFile": None,
+        "entFile": pathlib.Path(avaDir, "Inputs", "ENT", "entAlr.shp"),
+        "resFile": pathlib.Path(avaDir, "Inputs", "ENT", "entAlr.shp"),
     }
     variationDict = {"rho": np.asarray([200.0, 150.0]), "releaseScenario": ["relAlr"]}
 
     # call function to be tested
     simDict = com1DFA.prepareVarSimDict(standardCfg, inputSimFiles, variationDict)
+
     testCfg = configparser.ConfigParser()
     testCfg.optionxform = str
     testCfg["GENERAL"] = {
@@ -1755,6 +1953,8 @@ def test_prepareVarSimDict(tmp_path, caplog):
     }
     testCfg["INPUT"]["DEM"] = "avaAlr.tif"
     testCfg["INPUT"]["relThFile"] = ""
+    testCfg["INPUT"]["entrainment"] = str(pathlib.Path("ENT", "entAlr.shp"))
+    testCfg["INPUT"]["resistance"] = str(pathlib.Path("RES", "entAlr.shp"))
     testCfg["GENERAL"]["avalancheDir"] = str(avaDirTest)
 
     simHash = cfgUtils.cfgHash(testCfg)
@@ -1785,6 +1985,8 @@ def test_prepareVarSimDict(tmp_path, caplog):
         "entResInfo": {"flagEnt": "Yes", "flagRes": "Yes"},
         "demFile": avaDEM,
         "damFile": relPath,
+        "entFile": pathlib.Path(avaDir, "Inputs", "ENT", "entAlr.shp"),
+        "resFile": pathlib.Path(avaDir, "Inputs", "RES", "entAlr.shp"),
     }
     variationDict = {
         "rho": np.asarray([200.0, 150.0]),
@@ -1799,6 +2001,8 @@ def test_prepareVarSimDict(tmp_path, caplog):
         "entResInfo": {"flagEnt": "Yes", "flagRes": "Yes"},
         "demFile": avaDEM,
         "damFile": relPath,
+        "entFile": pathlib.Path(avaDir, "Inputs", "ENT", "entAlr.shp"),
+        "resFile": pathlib.Path(avaDir, "Inputs", "ENT", "entAlr.shp"),
     }
     testCfg2 = configparser.ConfigParser()
     testCfg2.optionxform = str
@@ -1839,10 +2043,12 @@ def test_prepareVarSimDict(tmp_path, caplog):
         "entThId": "0",
         "entThCi95": "None",
         "releaseScenario": "relAlr",
-        "DAM": str(pathlib.Path('DAM', relPath.name)),
+        "DAM": str(pathlib.Path("DAM", relPath.name)),
     }
     testCfg2["INPUT"]["DEM"] = "avaAlr.tif"
     testCfg2["INPUT"]["relThFile"] = ""
+    testCfg2["INPUT"]["entrainment"] = str(pathlib.Path("ENT", "entAlr.shp"))
+    testCfg2["INPUT"]["resistance"] = str(pathlib.Path("RES", "entAlr.shp"))
     testCfg2["GENERAL"]["avalancheDir"] = str(avaDirTest)
     simHash2 = cfgUtils.cfgHash(testCfg2)
     simName2 = "relAlr_" + simHash2 + "_C_L_entres_dfa"
@@ -1928,7 +2134,7 @@ def test_initializeSimulation(tmp_path):
     demOri = {"header": demHeader, "rasterData": demData}
 
     # setup release line, entrainment line
-    relFileTest = avaDir / 'REL' / "relAlr.shp"
+    relFileTest = avaDir / "REL" / "relAlr.shp"
     releaseLine = {
         "x": np.asarray([6.9, 8.5, 8.5, 6.9, 6.9]),
         "y": np.asarray([7.9, 7.9, 9.5, 9.5, 7.9]),
@@ -1941,7 +2147,7 @@ def test_initializeSimulation(tmp_path):
         "file": relFileTest,
     }
     entLine = {
-        "fileName": (avaDir / 'ENT' / "entAlr.shp"),
+        "fileName": (avaDir / "ENT" / "entAlr.shp"),
         "Name": ["testEnt"],
         "Start": np.asarray([0.0]),
         "thickness": [0.3, 0.3],
@@ -1977,7 +2183,10 @@ def test_initializeSimulation(tmp_path):
     #    print("reportAreaInfo", reportAreaInfo)
 
     pEnt = -10.0 * 2050.0 + 9.81 * 1.0
-    assert np.array_equal(particles["y"], np.asarray([6.25, 6.25, 6.25, 6.75, 7.25, 6.75, 6.75, 7.25, 7.25]))
+    assert np.array_equal(
+        particles["y"],
+        np.asarray([6.25, 6.25, 6.25, 6.75, 7.25, 6.75, 6.75, 7.25, 7.25]),
+    )
     assert np.sum(fields["pfv"]) == 0.0
     assert np.sum(fields["pft"]) != 0.0
     assert dem["header"]["xllcenter"] == 0.0
@@ -1986,10 +2195,17 @@ def test_initializeSimulation(tmp_path):
     assert dem["originalHeader"]["yllcenter"] == 2.0
     assert particles["nPart"] == 9
     assert np.array_equal(
-        particles["totalEnthalpy"], np.asarray([pEnt, pEnt, pEnt, pEnt, pEnt, pEnt, pEnt, pEnt, pEnt])
+        particles["totalEnthalpy"],
+        np.asarray([pEnt, pEnt, pEnt, pEnt, pEnt, pEnt, pEnt, pEnt, pEnt]),
+    )
+    assert np.array_equal(
+        particles["x"],
+        np.asarray([6.25, 6.75, 7.25, 6.25, 6.25, 6.75, 7.25, 6.75, 7.25]),
+    )
+    assert np.array_equal(
+        particles["m"],
+        np.asarray([50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0]),
     )
-    assert np.array_equal(particles["x"], np.asarray([6.25, 6.75, 7.25, 6.25, 6.25, 6.75, 7.25, 6.75, 7.25]))
-    assert np.array_equal(particles["m"], np.asarray([50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0]))
     assert particles["mTot"] == 450.0
     assert np.sum(particles["ux"]) == 0.0
     assert reportAreaInfo["Release area info"]["Projected Area [m2]"] == "4.00"
@@ -2004,7 +2220,7 @@ def test_initializeSimulation(tmp_path):
         "Start": np.asarray([0]),
         "Length": np.asarray([5]),
         "type": "Secondary release",
-        "fileName": (avaDir / 'SECREL' / "ec1.shp"),
+        "fileName": (avaDir / "SECREL" / "ec1.shp"),
         "Name": ["secRel1"],
         "thickness": [0.5],
         "thicknessSource": ["ini File"],
@@ -2040,7 +2256,6 @@ def test_initializeSimulation(tmp_path):
     assert reportAreaInfo["resistance"] == "No"
     assert np.sum(particles2["secondaryReleaseInfo"]["rasterData"]) == 4.5
 
-
     # test if dam is found
     # setup required input
     cfg = configparser.ConfigParser()
@@ -2074,7 +2289,7 @@ def test_initializeSimulation(tmp_path):
         "TIni": "-10.",
         "ResistanceModel": "cRes",
         "restitutionCoefficient": 1,
-        "nIterDam": 1
+        "nIterDam": 1,
     }
     releaseLine = {
         "x": np.asarray([6.9, 8.5, 8.5, 6.9, 6.9]),
@@ -2098,33 +2313,34 @@ def test_initializeSimulation(tmp_path):
         "muFile": None,
         "xiFile": None,
     }
-    inputSimLines['damLine'] = {
-        "fileName": [(avaDir / 'DAM' / "damLine.shp")],
-        'Name': [''],
-        'thickness': ['None'],
-        'slope': 60.0,
-        'Start': np.asarray([0.]),
-        'Length': np.asarray([2.]),
-        'x': np.asarray([5., 7.]),
-        'y': np.asarray([4., 6.]),
-        'z': np.asarray([0., 0.]),
-        'id': ['0'],
-        'ci95': ['None'],
-        'layerName': [None],
-        'nParts': [[0, 2]],
-        'nFeatures': 1,
-        'type': 'Dam'
+    inputSimLines["damLine"] = {
+        "fileName": [avaDir / "DAM" / "damLine.shp"],
+        "Name": [""],
+        "thickness": ["None"],
+        "slope": 60.0,
+        "Start": np.asarray([0.0]),
+        "Length": np.asarray([2.0]),
+        "x": np.asarray([5.0, 7.0]),
+        "y": np.asarray([4.0, 6.0]),
+        "z": np.asarray([0.0, 0.0]),
+        "id": ["0"],
+        "ci95": ["None"],
+        "layerName": [None],
+        "nParts": [[0, 2]],
+        "nFeatures": 1,
+        "type": "Dam",
     }
 
     particles3, fields3, dem3, reportAreaInfo3 = com1DFA.initializeSimulation(
         cfg, outDir, demOri, inputSimLines, logName
     )
 
-    print('dam', inputSimLines['damLine'])
+    print("dam", inputSimLines["damLine"])
+
+    assert reportAreaInfo3["dam"] == "Yes"
+    assert "xCrown" in inputSimLines["damLine"]
+    assert "height" in inputSimLines["damLine"]
 
-    assert reportAreaInfo3['dam'] == 'Yes'
-    assert 'xCrown' in inputSimLines['damLine']
-    assert 'height' in inputSimLines['damLine']
 
 def test_runCom1DFA(tmp_path, caplog):
     """Check that runCom1DFA produces the good outputs"""
@@ -2195,6 +2411,7 @@ def test_runCom1DFA(tmp_path, caplog):
         "nExitedParticles",
         "dmDet",
         "massDetrained",
+        "tPlot",
     ]
 
     # read one particles dictionary
@@ -2313,7 +2530,10 @@ def test_fetchRelVolume(tmp_path):
     relThField1 = pathlib.Path(avaDir, "Inputs", "RELTH", "relThField1.asc")
     fU.makeADir(pathlib.Path(avaDir, "Inputs", "RELTH"))
     IOf.writeResultToRaster(
-        relThF["header"], relThF["rasterData"], relThField1.parent / relThField1.stem, flip=False
+        relThF["header"],
+        relThF["rasterData"],
+        relThField1.parent / relThField1.stem,
+        flip=False,
     )
 
     cfg = {}
@@ -2329,7 +2549,12 @@ def test_fetchRelVolume(tmp_path):
         "relTh1": 4.0,
         "secRelArea": False,
     }
-    cfg["INPUT"] = {"relThFile": "", "relThId": "0|1", "relThThickness": "2.|4.", "thFromIni": ""}
+    cfg["INPUT"] = {
+        "relThFile": "",
+        "relThId": "0|1",
+        "relThThickness": "2.|4.",
+        "thFromIni": "",
+    }
 
     relVolume = com1DFA.fetchRelVolume(rel1, cfg, demPath, None)
 
diff --git a/avaframe/tests/test_com1DFATools.py b/avaframe/tests/test_com1DFATools.py
index dc2b9f993..108535b79 100644
--- a/avaframe/tests/test_com1DFATools.py
+++ b/avaframe/tests/test_com1DFATools.py
@@ -1,5 +1,5 @@
 """
-    Pytest for module com1DFATools
+Pytest for module com1DFATools
 """
 
 #  Load modules
@@ -8,61 +8,145 @@
 import configparser
 import pathlib
 import shutil
+import numpy as np
+import rasterio
 
 from avaframe.com1DFA import com1DFATools
+from avaframe.in2Trans.rasterUtils import transformFromASCHeader
+import avaframe.in3Utils.fileHandlerUtils as fU
 
 
 def test_getPartInitMethod(capfd):
     cfg = configparser.ConfigParser()
     csz = 1
     relThForPart = 2
-    cfg['GENERAL'] = {'rho': '3', 'massPerPart': '10', 'deltaTh': '0.25', 'sphKernelRadius': '1', 'nPPK0': '5',
-                      'aPPK': '-1', 'sphKR0': '5', 'massPerParticleDeterminationMethod': 'MPPDIR'}
-    massPerPart, nPPK = com1DFATools.getPartInitMethod(cfg['GENERAL'], csz, relThForPart)
+    cfg["GENERAL"] = {
+        "rho": "3",
+        "massPerPart": "10",
+        "deltaTh": "0.25",
+        "sphKernelRadius": "1",
+        "nPPK0": "5",
+        "aPPK": "-1",
+        "sphKR0": "5",
+        "massPerParticleDeterminationMethod": "MPPDIR",
+    }
+    massPerPart, nPPK = com1DFATools.getPartInitMethod(cfg["GENERAL"], csz, relThForPart)
     assert massPerPart == 10
     assert nPPK == 0
 
-    cfg['GENERAL']['massPerParticleDeterminationMethod'] = 'MPPDH'
-    massPerPart, nPPK = com1DFATools.getPartInitMethod(cfg['GENERAL'], csz, relThForPart)
+    cfg["GENERAL"]["massPerParticleDeterminationMethod"] = "MPPDH"
+    massPerPart, nPPK = com1DFATools.getPartInitMethod(cfg["GENERAL"], csz, relThForPart)
     assert massPerPart == 0.75
     assert nPPK == 0
 
-    cfg['GENERAL']['massPerParticleDeterminationMethod'] = 'MPPKR'
-    massPerPart, nPPK = com1DFATools.getPartInitMethod(cfg['GENERAL'], csz, relThForPart)
+    cfg["GENERAL"]["massPerParticleDeterminationMethod"] = "MPPKR"
+    massPerPart, nPPK = com1DFATools.getPartInitMethod(cfg["GENERAL"], csz, relThForPart)
 
     assert massPerPart == pytest.approx(math.pi * 6 / 25, abs=1e-6)
     assert nPPK == 25
 
 
 def test_createSimDictFromCfgs(tmp_path):
-    """ test creating a simDict from multiple cfg files """
+    """test creating a simDict from multiple cfg files"""
 
     dirPath = pathlib.Path(__file__).parents[0]
-    testPath = dirPath / 'data' / 'com1DFAConfigsTest'
-    inputDir = dirPath / 'data' / 'testCom1DFA2'
-    avaDir = pathlib.Path(tmp_path, 'testCom1DFA')
+    testPath = dirPath / "data" / "com1DFAConfigsTest"
+    inputDir = dirPath / "data" / "testCom1DFA2"
+    avaDir = pathlib.Path(tmp_path, "testCom1DFA")
     shutil.copytree(inputDir, avaDir)
     cfgMain = configparser.ConfigParser()
-    cfgMain['MAIN'] = {'avalancheDir': avaDir}
+    cfgMain["MAIN"] = {"avalancheDir": avaDir}
 
     with pytest.raises(FileNotFoundError) as e:
         # call function to be tested
         simDict, inputSimFiles, simDFExisting, outDir = com1DFATools.createSimDictFromCfgs(cfgMain, testPath)
     assert ("No configuration file found to create simulation runs in:") in str(e.value)
 
-
-    testPath = dirPath / 'data' / 'com1DFAConfigs'
-    avaDir = pathlib.Path(tmp_path, 'testCom1DFA2')
+    testPath = dirPath / "data" / "com1DFAConfigs"
+    avaDir = pathlib.Path(tmp_path, "testCom1DFA2")
     shutil.copytree(inputDir, avaDir)
-    cfgMain['MAIN'] = {'avalancheDir': avaDir}
+    cfgMain["MAIN"] = {"avalancheDir": avaDir}
     # call function to be tested
     simDict, inputSimFiles, simDFExisting, outDir = com1DFATools.createSimDictFromCfgs(cfgMain, testPath)
 
     # for sim in simDict.keys():
-#        print('simDict: ', sim, simDict[sim])
-#    print('inputSimFiles', inputSimFiles)
-#    print('simDFExisting', simDFExisting)
-#    print('outDir', outDir)
+    #        print('simDict: ', sim, simDict[sim])
+    #    print('inputSimFiles', inputSimFiles)
+    #    print('simDFExisting', simDFExisting)
+    #    print('outDir', outDir)
 
     assert simDFExisting == None
     assert len(simDict) == 16
+
+
+def test_updateResCoeffFields(tmp_path):
+    """test updating fields of cResH and detK for resistance"""
+
+    outDir = pathlib.Path(tmp_path, "Outputs", "com1DFA", "reports")
+    fU.makeADir(outDir)
+    cfg = configparser.ConfigParser()
+    cfg["GENERAL"] = {
+        "forestVMin": 6.0,
+        "forestThMin": 0.6,
+        "forestVMax": 40.0,
+        "forestThMax": 10.0,
+        "avalancheDir": pathlib.Path(tmp_path),
+    }
+
+    cResRasterOrig = np.ones((4, 10))
+    detRasterOrig = np.ones((4, 10)) + 2
+    FV = np.zeros((4, 10))
+    FT = np.zeros((4, 10))
+
+    fields = {
+        "cResRasterOrig": cResRasterOrig,
+        "detRasterOrig": detRasterOrig,
+        "FV": FV,
+        "FT": FT,
+    }
+    dem = {
+        "originalHeader": {
+            "ncols": 10,
+            "nrows": 4,
+            "xllcenter": 1000.00,
+            "yllcenter": -5000.00,
+            "cellsize": 5,
+            "nodata_value": -9999.00,
+            "driver": "AAIGrid",
+        }
+    }
+    transform = transformFromASCHeader(dem["originalHeader"])
+    dem["originalHeader"]["transform"] = transform
+    dem["originalHeader"]["crs"] = rasterio.crs.CRS()
+
+    fields = com1DFATools.updateResCoeffFields(fields, cfg["GENERAL"], 0.0, dem)
+
+    assert fields["cResRaster"].any() == False
+
+    FV[0, 4] = 7
+    FT[0, 4] = 0.7
+    FT[2, 2] = 0.7
+    FV[1, 3] = 41
+    FT[1, 1] = 11
+    fields = {
+        "cResRasterOrig": cResRasterOrig,
+        "detRasterOrig": detRasterOrig,
+        "FV": FV,
+        "FT": FT,
+    }
+
+    fields = com1DFATools.updateResCoeffFields(fields, cfg["GENERAL"], 0.0, dem)
+
+    assert fields["cResRaster"].any() == True
+    assert fields["cResRaster"][0, 4] == 1
+    assert fields["detRaster"][0, 4] == 0
+    assert fields["cResRaster"][2, 2] == 0
+    assert fields["detRaster"][2, 2] == 3
+    assert fields["cResRaster"][1, 3] == 0
+    assert fields["detRaster"][1, 3] == 0
+    assert fields["cResRaster"][1, 1] == 0
+    assert fields["detRaster"][1, 1] == 0
+    assert fields["cResRaster"][0, 0] == 0
+    assert fields["detRaster"][0, 0] == 3
+    assert fields["cResRaster"][2, 3] == 0
+    assert fields["detRaster"][2, 3] == 3
diff --git a/docs/moduleCom1DFA.rst b/docs/moduleCom1DFA.rst
index 3a8d6de29..26a87f396 100644
--- a/docs/moduleCom1DFA.rst
+++ b/docs/moduleCom1DFA.rst
@@ -279,6 +279,7 @@ The result types that can be chosen to be exported are (all correspond to fields
 * Vx, Vy, Vz - velocity x-, y- and z-component
 * TA - travel angle
 * dmDet - detrained mass
+* FTDet - thickness of detrained mass computed based on dmDet / (rho * area of cell)
 * particles (:ref:`com1DFAAlgorithm:Particle properties`)
 
 Have a look at the designated subsection Output in ``com1DFA/com1DFACfg.ini``.