diff --git a/_quarto.yml b/_quarto.yml
index 606e6d9da7c95e1620eaa86fee23020f8588b220..0fe852244b687e738c476241e6d374a50c7c87b2 100644
--- a/_quarto.yml
+++ b/_quarto.yml
@@ -52,7 +52,7 @@ website:
- "exercises/user-experience.qmd"
# - "exercises/testing.qmd"
# - "exercises/git2.qmd"
- - "exercises/parallelism.qmd"
+ - "exercises/parallelism/parallelism.qmd"
# - "exercises/hardware.qmd"
# - "exercises/file_and_data_systems.qmd"
# - "exercises/memory_hierarchies.qmd"
diff --git a/exercises/parallelism/main.c b/exercises/parallelism/main.c
new file mode 100644
index 0000000000000000000000000000000000000000..f4a474945b39d363a17dcddfd727dd75936de571
--- /dev/null
+++ b/exercises/parallelism/main.c
@@ -0,0 +1,297 @@
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+
+#ifdef _OPENMP
+#include <omp.h>
+#else
+#ifdef WRITE_DECOMP
+#error You cannot write the decomposition without OpenMP.
+#endif
+#endif
+
+//compile using (for example):
+//gcc -lm -fopenmp -g main.c
+
+//definition of constants
+//some values are given in meter
+//some in seconds
+//
+#define x_0 10000
+#define delta_x 1
+#define coeff_sech 0.01
+//depth of sea at deepest point
+#define b_0 1000
+//gravitational acceleration
+#define g 10
+#define delta_t 0.001
+
+
+double *h;
+double *h1;
+double *h2;
+double *u;
+double *u1;
+double *u2;
+
+int *point_to_thread;
+
+long count_points () {
+ return (int)(2*x_0/delta_x);
+}
+
+double get_x(long i){
+ return i*delta_x - x_0;
+}
+
+double waterdepth(long i){
+ double x = get_x(i);
+ return b_0*(x*x/(x_0*x_0)-1);
+}
+
+void initialise_h(double *h){
+ double x, tmp;
+ long i;
+ for(i=0; i < count_points(); i++){
+ x=get_x(i);
+ tmp = 2/(exp(coeff_sech*x)+exp(-coeff_sech*x));
+ if (tmp > 1e-10)
+ h[i] = tmp*tmp;
+ else
+ h[i] = 0;
+ }
+}
+
+void initialise_u(double *u){
+ long i;
+ for(i=0;i<count_points();i++){
+ u[i]=0;
+ }
+}
+
+
+
+void print_values(FILE *data, double *h, double *u){
+ long i;
+ fprintf(data, "#x\th\tu\tb\n");
+ for(i=0;i<count_points();i++){
+ fprintf(data, "%g\t%g\t%g\t%g\n", get_x(i), h[i], u[i], waterdepth(i));
+ }
+}
+
+double lu(double *h, double *u, long i){
+ return u[i]*(u[i+1]-u[i-1])+g*(h[i+1]-h[i-1]);
+}
+
+double lh(double *h, double *u, long i){
+ return u[i]*(h[i+1]-h[i-1])-u[i]*(waterdepth(i+1)-waterdepth(i-1))+(h[i]-waterdepth(i))*(u[i+1]-u[i-1]);
+}
+
+
+//sets h1,u1 using values (startconditions) from u2/h2
+void euler(double *h1, double *h2, double *u1, double *u2){
+ long i;
+ long npoints = count_points();
+ u1[0]=0; //boundary value
+ h1[0]=0; //boundary value
+ u1[npoints-1]=0;
+ h1[npoints-1]=0;
+ for(i=1;i<npoints-1;i++){
+ u1[i]=u2[i]-(delta_t/(2.*delta_x))*lu(h2,u2,i);
+ h1[i]=h2[i]-(delta_t/(2.*delta_x))*lh(h2,u2,i);
+ }
+}
+
+long findmax(double *h, long oldmaxi,long outputmth){
+ long i, maxi=oldmaxi-10;
+ for(i=oldmaxi-10;(i< oldmaxi+(int)(outputmth/delta_x)+30) && (i<count_points()) ;i++){
+ if(h[i]>h[maxi]){
+ maxi=i;
+ }
+ }
+ return maxi;
+}
+
+double maxval(double *h) {
+ long npoints = count_points();
+ long i;
+ double hmax = -10;
+// HOMEWORK: Parallelise this loop!
+ for(i = 0; i < npoints; i++){
+ if (h[i] > hmax)
+ hmax = h[i];
+ }
+ return hmax;
+}
+
+
+double ng_derivative(double xi, double xi_1, double deltax){
+ return (xi-xi_1)/deltax;
+}
+
+
+void siginthandler(){
+ print_values(stdout,h,u);
+ exit(0);
+}
+
+#ifdef _OPENMP
+void save_decomp(long i) {
+ point_to_thread[i] = omp_get_thread_num();
+}
+
+void print_decomp() {
+ int i,j;
+ int npoints = count_points();
+ for (i=0; i<omp_get_num_threads(); i++) {
+ #pragma omp barrier
+ if (i != omp_get_thread_num())
+ continue;
+ printf("%d\t", i);
+ int first=-1;
+ for (j=0; j<npoints+1; j++) {
+ if (j < npoints && point_to_thread[j] == i) {
+ if (first == -1)
+ first = j;
+ } else {
+ if (first == -1)
+ continue;
+ if (first == j-1)
+ printf("%d ", j);
+ else {
+ printf("%d-%d ", first, j-1);
+ }
+ first = -1;
+ }
+ }
+ printf("\n");
+ }
+}
+#endif
+
+
+void swap(double **var, double **var1, double **var2) {
+ double *temp;
+
+ temp = *var2;
+ *var2 = *var1;
+ *var1 = *var;
+ *var = temp;
+}
+
+void set_boundaries(double *var) {
+ long npoints = count_points();
+ var[0] = 0;
+ var[npoints-1] = 0;
+}
+
+int main(int argc, char* argv){
+ long i;
+ long npoints = count_points();
+ h = malloc(sizeof(double)*npoints);
+ h1 = malloc(sizeof(double)*npoints);
+ h2 = malloc(sizeof(double)*npoints);
+
+ u = malloc(sizeof(double)*npoints);
+ u1 = malloc(sizeof(double)*npoints);
+ u2 = malloc(sizeof(double)*npoints);
+#ifdef WRITE_DECOMP
+ point_to_thread = malloc(sizeof(int)*count_points());
+#endif
+
+ if(h==NULL || h1==NULL || h2==NULL || u==NULL || u1==NULL || u2==NULL){
+ fprintf(stderr, "Insufficient memory");
+ exit(23);
+ }
+#if 0
+ signal(SIGINT, siginthandler);
+#endif
+
+ long maxi=((npoints)/2)+5;
+ long oldmaxi;
+ long oldmaxi_output=(npoints)/2;
+ int maximth=20;
+ double t=0;
+ long n=0;
+ int outputmth=1000; //write every mth set of values to stderr
+ int dataqth=200; //write every qth set of values to output#
+ char file[60]; //char pointer for filename
+ FILE* out;
+ double *temp;
+ int first = 1;
+
+ initialise_u(u2);
+ initialise_h(h2);
+ euler(h1,h2,u1,u2);
+#if 0
+ out = fopen("t-1.dat","w");
+ print_values(out,h1,u1);
+ out = fopen("t-2.dat","w");
+ print_values(out,h2,u2);
+#endif
+ time_t t_before = time(NULL);
+ fprintf(stderr, "#time h_max\n");
+
+ //while ( t < 150 ){
+ while ( t < 50 ){
+ if (!first) {
+ swap(&h, &h1, &h2);
+ swap(&u, &u1, &u2);
+ } else
+ first = 0;
+
+ set_boundaries(u);
+ set_boundaries(h);
+
+ #pragma omp parallel for // schedule(dynamic, 10)
+ for(i=0;i<npoints;i++){
+#ifdef WRITE_DECOMP
+ save_decomp(i);
+#endif
+ if (i == 0 || i == npoints-1)
+ continue;
+ u[i]=u2[i]-(delta_t/delta_x)*lu(h1,u1,i);
+ h[i]=h2[i]-(delta_t/delta_x)*lh(h1,u1,i);
+ }
+
+ t+=delta_t;
+ if((n%maximth)==0){
+
+ oldmaxi=maxi;
+ maxi=findmax(h,oldmaxi,maximth);
+ if((n%outputmth)==0){
+ double maxh = maxval(h);
+ fprintf(stderr, "%g\t%g\n", t, maxh);
+
+#if 0
+ fprintf(stderr, "%g \t",t);
+ fprintf(stderr, "%g\t%g \t%g\n",get_x(maxi),h[maxi],ng_derivative(get_x(maxi),get_x(oldmaxi_output),outputmth*delta_t));
+ oldmaxi_output=maxi;
+#endif
+ }
+ }
+
+
+#if 0
+ if((n%dataqth)==0){
+ sprintf(file,"output%04ld.dat",n/dataqth);
+ out = fopen(file,"w");
+ print_values(out,h,u);
+ }
+#endif
+ n++;
+ }
+ time_t t_after = time(NULL);
+ double duration = difftime(t_after, t_before);
+ printf("Run took %g seconds.\n", duration);
+
+ //print_values(stdout,h,u);
+#ifdef WRITE_DECOMP
+#pragma omp parallel
+ print_decomp();
+#endif
+
+ return 0;
+}
+
diff --git a/exercises/parallelism.qmd b/exercises/parallelism/parallelism.qmd
similarity index 100%
rename from exercises/parallelism.qmd
rename to exercises/parallelism/parallelism.qmd