diff --git a/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp b/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp
index 4ef63c614c38cd4f6ea15179fff4aadcfb8b3e6a..bc414530f64193223e63e673412c67bfc6aef745 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp
+++ b/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp
@@ -38,7 +38,7 @@ void CudaMemoryManager::cudaAllocCoord(int lev)
 	checkCudaErrors( cudaMallocHost((void**) &(parameter->getParH(lev)->coordX_SP      ), parameter->getParH(lev)->mem_size_real_SP  ));
 	checkCudaErrors( cudaMallocHost((void**) &(parameter->getParH(lev)->coordY_SP      ), parameter->getParH(lev)->mem_size_real_SP  ));
 	checkCudaErrors( cudaMallocHost((void**) &(parameter->getParH(lev)->coordZ_SP      ), parameter->getParH(lev)->mem_size_real_SP  ));
-	//Device (spinning ship)
+	//Device (spinning ship + upsala)
 	checkCudaErrors( cudaMalloc((void**) &(parameter->getParD(lev)->coordX_SP      ), parameter->getParH(lev)->mem_size_real_SP  ));
 	checkCudaErrors( cudaMalloc((void**) &(parameter->getParD(lev)->coordY_SP      ), parameter->getParH(lev)->mem_size_real_SP  ));
 	checkCudaErrors( cudaMalloc((void**) &(parameter->getParD(lev)->coordZ_SP      ), parameter->getParH(lev)->mem_size_real_SP  ));
@@ -59,6 +59,36 @@ void CudaMemoryManager::cudaFreeCoord(int lev)
 	checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->coordY_SP   ));
 	checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->coordZ_SP   ));
 }
+void CudaMemoryManager::cudaAllocBodyForce(int lev) 
+{
+    //Host
+	checkCudaErrors( cudaMallocHost((void**) &(parameter->getParH(lev)->forceX_SP      ), parameter->getParH(lev)->mem_size_real_SP  ));
+	checkCudaErrors( cudaMallocHost((void**) &(parameter->getParH(lev)->forceY_SP      ), parameter->getParH(lev)->mem_size_real_SP  ));
+	checkCudaErrors( cudaMallocHost((void**) &(parameter->getParH(lev)->forceZ_SP      ), parameter->getParH(lev)->mem_size_real_SP  ));
+	//Device (spinning ship)
+	checkCudaErrors( cudaMalloc((void**) &(parameter->getParD(lev)->forceX_SP      ), parameter->getParH(lev)->mem_size_real_SP  ));
+	checkCudaErrors( cudaMalloc((void**) &(parameter->getParD(lev)->forceY_SP      ), parameter->getParH(lev)->mem_size_real_SP  ));
+	checkCudaErrors( cudaMalloc((void**) &(parameter->getParD(lev)->forceZ_SP      ), parameter->getParH(lev)->mem_size_real_SP  ));
+	//////////////////////////////////////////////////////////////////////////
+	double tmp = 3. * (double)parameter->getParH(lev)->mem_size_real_SP;
+	setMemsizeGPU(tmp, false);
+
+}
+void CudaMemoryManager::cudaCopyBodyForce(int lev) 
+{
+   	//copy host to device
+	checkCudaErrors( cudaMemcpy(parameter->getParD(lev)->forceX_SP,  parameter->getParH(lev)->forceX_SP,  parameter->getParH(lev)->mem_size_real_SP     , cudaMemcpyHostToDevice));
+	checkCudaErrors( cudaMemcpy(parameter->getParD(lev)->forceY_SP,  parameter->getParH(lev)->forceY_SP,  parameter->getParH(lev)->mem_size_real_SP     , cudaMemcpyHostToDevice));
+	checkCudaErrors( cudaMemcpy(parameter->getParD(lev)->forceZ_SP,  parameter->getParH(lev)->forceZ_SP,  parameter->getParH(lev)->mem_size_real_SP     , cudaMemcpyHostToDevice));
+
+}
+void CudaMemoryManager::cudaFreeBodyForce(int lev) 
+{
+   	checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->forceX_SP   ));
+	checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->forceY_SP   ));
+	checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->forceZ_SP   ));
+
+}
 //print
 void CudaMemoryManager::cudaCopyDataToHost(int lev)
 {
diff --git a/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.h b/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.h
index 315c40e5d9dd0784a9753d565e9159449c396a33..f96b93ba39ed9c37efa7922373d91d30f0cb5eb8 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.h
+++ b/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.h
@@ -37,7 +37,11 @@ public:
 	void cudaCopyCoord(int lev);
 	void cudaFreeCoord(int lev);
 
-	void cudaCopyDataToHost(int lev);
+	void cudaAllocBodyForce(int lev);
+    void cudaCopyBodyForce(int lev);
+    void cudaFreeBodyForce(int lev);
+
+    void cudaCopyDataToHost(int lev);
 
 	void cudaAllocSP(int lev);
 	void cudaCopySP(int lev);
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
index d0a99c9f718145fb6764b3ef044c47fcc389f4d0..89a50c7272d42516fe23824f7686284773a6fbe0 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
@@ -3666,7 +3666,12 @@ void Parameter::setSimulatePorousMedia(bool simulatePorousMedia)
 
 void Parameter::setIsF3(bool isF3)
 {
-	this->isF3 = isF3;
+	this->isF3 = isF3; 
+}
+
+void Parameter::setIsBodyForce(bool isBodyForce) 
+{
+	this->isBodyForce = isBodyForce;
 }
 
 void Parameter::setGridX(std::vector<int> GridX)
@@ -5017,7 +5022,12 @@ bool Parameter::getSimulatePorousMedia()
 
 bool Parameter::getIsF3()
 {
-	return this->isF3;
+	return this->isF3; 
+}
+
+bool Parameter::getIsBodyForce() 
+{ 
+	return this->isBodyForce; 
 }
 
 bool Parameter::getIsGeometryValues()
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
index c7696eac31203a8bf26724141071ad84ba27bcda..e806246676ea126507ff2a0c0d44d875c6440a4c 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
@@ -80,6 +80,9 @@ struct ParameterStruct{
 	//unsigned int *coordX_SP, *coordY_SP, *coordZ_SP;
 	real *coordX_SP, *coordY_SP, *coordZ_SP;
 
+	//body forces////////////
+	real *forceX_SP, *forceY_SP, *forceZ_SP;
+
 	//vel parab///////////////
 	real *vParab;
 
@@ -723,6 +726,7 @@ public:
 	void setUseInitNeq(bool useInitNeq);
 	void setSimulatePorousMedia(bool simulatePorousMedia);
 	void setIsF3(bool isF3);
+    void setIsBodyForce(bool isBodyForce);
 	void setclockCycleForMP(real clockCycleForMP);
 	void setDevices(std::vector<uint> devices);
 	void setGridX(std::vector<int> GridX);
@@ -972,6 +976,7 @@ public:
 	bool getUseInitNeq();
 	bool getSimulatePorousMedia();
 	bool getIsF3();
+    bool getIsBodyForce();
 	double getMemsizeGPU();
 	//1D domain decomposition
 	std::vector<std::string> getPossNeighborFiles(std::string sor);
@@ -1033,6 +1038,7 @@ private:
 	bool calcDragLift, calcCp;
 	bool writeVeloASCII;
 	bool calcPlaneConc;
+    bool isBodyForce;
 	int diffMod;
 	int coarse, fine, maxlevel;
 	int factor_gridNZ;