diff --git a/lib/calc/xnmedian.c b/lib/calc/xnmedian.c
index 3dab302e90e..e10efbce1ac 100644
--- a/lib/calc/xnmedian.c
+++ b/lib/calc/xnmedian.c
@@ -44,8 +44,7 @@ static int dcmp(const void *aa, const void *bb)
 
 int f_nmedian(int argc, const int *argt, void **args)
 {
-    static void *array;
-    static int alloc;
+    void *array;
     int size = argc * Rast_cell_size(argt[0]);
     int i, j;
 
@@ -56,10 +55,7 @@ int f_nmedian(int argc, const int *argt, void **args)
         if (argt[i] != argt[0])
             return E_ARG_TYPE;
 
-    if (size > alloc) {
-        alloc = size;
-        array = G_realloc(array, size);
-    }
+    array = G_malloc(size);
 
     switch (argt[0]) {
     case CELL_TYPE: {
diff --git a/raster/r.mapcalc/CMakeLists.txt b/raster/r.mapcalc/CMakeLists.txt
index ab9155ce0bb..e851137e2b5 100644
--- a/raster/r.mapcalc/CMakeLists.txt
+++ b/raster/r.mapcalc/CMakeLists.txt
@@ -30,7 +30,8 @@ build_program(
   OPTIONAL_DEPENDS
   Readline::Readline
   Readline::History
-  Threads::Threads)
+  Threads::Threads
+  OPENMP)
 
 build_program(
   NAME
diff --git a/raster/r.mapcalc/Makefile b/raster/r.mapcalc/Makefile
index 7b8bbf90d16..d8ed91eeb1c 100644
--- a/raster/r.mapcalc/Makefile
+++ b/raster/r.mapcalc/Makefile
@@ -12,8 +12,10 @@ r3_mapcalc_OBJS := $(filter-out map.o xcoor.o xres.o, $(AUTO_OBJS))
 include $(MODULE_TOPDIR)/include/Make/Multi.make
 
 EXTRA_CFLAGS = $(READLINEINCPATH) $(PTHREADINCPATH)
-LIBES2 = $(CALCLIB) $(GISLIB) $(RASTERLIB) $(BTREELIB) $(READLINELIBPATH) $(READLINELIB) $(HISTORYLIB) $(PTHREADLIBPATH) $(PTHREADLIB)
-LIBES3 = $(CALCLIB) $(RASTER3DLIB) $(GISLIB) $(RASTERLIB) $(BTREELIB) $(READLINELIBPATH) $(READLINELIB) $(HISTORYLIB) $(PTHREADLIBPATH) $(PTHREADLIB)
+LIBES2 = $(CALCLIB) $(GISLIB) $(RASTERLIB) $(BTREELIB) $(READLINELIBPATH) $(READLINELIB) $(HISTORYLIB) $(PTHREADLIBPATH) $(PTHREADLIB) $(OPENMP_LIBPATH) $(OPENMP_LIB)
+LIBES3 = $(CALCLIB) $(RASTER3DLIB) $(GISLIB) $(RASTERLIB) $(BTREELIB) $(READLINELIBPATH) $(READLINELIB) $(HISTORYLIB) $(PTHREADLIBPATH) $(PTHREADLIB) $(OPENMP_LIBPATH) $(OPENMP_LIB)
+EXTRA_CFLAGS = $(OPENMP_CFLAGS)
+EXTRA_INC = $(OPENMP_INCPATH)
 
 default: multi
 
diff --git a/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py b/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py
new file mode 100644
index 00000000000..5cd747f28a4
--- /dev/null
+++ b/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py
@@ -0,0 +1,104 @@
+"""Benchmarking of r.mapcalc
+raster (2D)
+
+@author Chung-Yuan Liang, 2025
+"""
+
+from grass.exceptions import CalledModuleError
+from grass.pygrass.modules import Module
+
+import grass.benchmark as bm
+
+
+def main():
+    results = []
+    metrics = ["time", "speedup", "efficiency"]
+    mapsizes = [10e6, 50e6, 100e6]
+
+    # run benchmarks
+
+    for mapsize in mapsizes:
+        benchmark(
+            size=int(mapsize**0.5),
+            step=0,
+            results=results,
+        )
+
+    # plot results
+
+    for metric in metrics:
+        bm.nprocs_plot(
+            results,
+            filename=f"r_mapcalc_{metric}.svg",
+            title=f"r.mapcalc {metric}",
+            metric=metric,
+        )
+
+
+def benchmark(size, step, results):
+    map1 = "benchmark_r_mapcalc_map1"
+    map2 = "benchmark_r_mapcalc_map2"
+    output = "benchmark_r_mapcalc"
+
+    generate_map(rows=size, cols=size, fname=map1)
+    generate_map(rows=size, cols=size, fname=map2)
+    module = Module(
+        "r.mapcalc",
+        expression=f"{output}=({map1} + {map2})",
+        overwrite=True,
+    )
+
+    results.append(
+        bm.benchmark_nprocs(
+            module,
+            label=f"r.mapcalc_simple_{int((size * size) / 1e6)}M",
+            max_nprocs=12,
+            repeat=5,
+        )
+    )
+
+    module = Module(
+        "r.mapcalc",
+        expression=f"{output}=({map1} + {map2} - 2*{map2} + {map1}*{map2} - {map1}/{map2})/2",
+        overwrite=True,
+    )
+    results.append(
+        bm.benchmark_nprocs(
+            module,
+            label=f"r.mapcalc_complex_{int((size * size) / 1e6)}M",
+            max_nprocs=12,
+            repeat=5,
+        )
+    )
+
+    module = Module(
+        "r.mapcalc",
+        expression=f"{output}= if (({map1}[5, 5] + {map2}[-5, -5]) > 0, 1.6, 0.1)",
+        overwrite=True,
+    )
+    results.append(
+        bm.benchmark_nprocs(
+            module,
+            label=f"r.mapcalc_neighbor_{int((size * size) / 1e6)}M",
+            max_nprocs=12,
+            repeat=5,
+        )
+    )
+    Module(
+        "g.remove", quiet=True, flags="f", type="raster", pattern="benchmark_r_mapcalc*"
+    )
+
+
+def generate_map(rows, cols, fname):
+    Module("g.region", flags="p", rows=rows, cols=cols, res=1)
+    # Generate using r.random.surface if r.surf.fractal fails
+    try:
+        print("Generating reference map using r.surf.fractal...")
+        Module("r.surf.fractal", output=fname, overwrite=True)
+    except CalledModuleError:
+        print("r.surf.fractal fails, using r.random.surface instead...")
+        Module("r.random.surface", output=fname, overwrite=True)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/raster/r.mapcalc/evaluate.c b/raster/r.mapcalc/evaluate.c
index e5c2dcdd9d7..a0bac5f9cc1 100644
--- a/raster/r.mapcalc/evaluate.c
+++ b/raster/r.mapcalc/evaluate.c
@@ -1,3 +1,7 @@
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
+
 #include <stdlib.h>
 #include <unistd.h>
 #include <string.h>
@@ -12,8 +16,9 @@
 
 /****************************************************************************/
 
-int current_depth, current_row;
-int depths, rows, columns;
+int current_depth;
+int *current_row;
+int depths, rows;
 
 /* Local variables for map management */
 static expression **map_list = NULL;
@@ -69,14 +74,50 @@ void extract_maps(expression *e)
 
 static void allocate_buf(expression *e)
 {
-    e->buf = G_malloc(columns * Rast_cell_size(e->res_type));
+
+    int threads = 1;
+#if defined(_OPENMP)
+    threads = omp_get_max_threads();
+#endif
+
+    e->buf = (void **)G_malloc(sizeof(void *) * threads);
+    for (int t = 0; t < threads; t++)
+        e->buf[t] = G_malloc(columns * Rast_cell_size(e->res_type));
 }
 
-static void set_buf(expression *e, void *buf)
+static void set_buf(expression *e, void **buf)
 {
     e->buf = buf;
 }
 
+static void free_buf(expression *e)
+{
+    int threads = 1;
+#if defined(_OPENMP)
+    threads = omp_get_max_threads();
+#endif
+
+    for (int t = 0; t < threads; t++) {
+        G_free(e->buf[t]);
+        e->buf[t] = NULL;
+    }
+    G_free(e->buf);
+    e->buf = NULL;
+}
+
+static void free_argv(expression *e)
+{
+    int i;
+
+    for (i = 1; i <= e->data.func.argc; i++) {
+        free_buf(e->data.func.args[i]);
+        e->data.func.args[i]->buf = NULL;
+    }
+
+    G_free(e->data.func.argv);
+    e->data.func.argv = NULL;
+}
+
 /****************************************************************************/
 
 static void initialize_constant(expression *e)
@@ -93,8 +134,15 @@ static void initialize_map(expression *e)
 {
     allocate_buf(e);
 
-    e->data.map.idx = open_map(e->data.map.name, e->data.map.mod,
-                               e->data.map.row, e->data.map.col);
+    int threads = 1;
+#if defined(_OPENMP)
+    threads = omp_get_max_threads();
+#endif
+    e->data.map.idx = G_malloc(threads * sizeof(int));
+    for (int t = 0; t < threads; t++) {
+        e->data.map.idx[t] = open_map(e->data.map.name, e->data.map.mod,
+                                      e->data.map.row, e->data.map.col);
+    }
 }
 
 static void initialize_function(expression *e)
@@ -102,8 +150,7 @@ static void initialize_function(expression *e)
     int i;
 
     allocate_buf(e);
-
-    e->data.func.argv = G_malloc((e->data.func.argc + 1) * sizeof(void *));
+    e->data.func.argv = G_malloc((e->data.func.argc + 1) * sizeof(void **));
     e->data.func.argv[0] = e->buf;
 
     for (i = 1; i <= e->data.func.argc; i++) {
@@ -163,9 +210,14 @@ static void end_evaluate(struct expression *e)
 
 static void evaluate_constant(expression *e)
 {
-    int *ibuf = e->buf;
-    float *fbuf = e->buf;
-    double *dbuf = e->buf;
+    int tid = 0;
+#if defined(_OPENMP)
+    tid = omp_get_thread_num();
+#endif
+
+    int *ibuf = e->buf[tid];
+    float *fbuf = e->buf[tid];
+    double *dbuf = e->buf[tid];
     int i;
 
     switch (e->res_type) {
@@ -195,15 +247,24 @@ static void evaluate_variable(expression *e UNUSED)
 
 static void evaluate_map(expression *e)
 {
-    get_map_row(
-        e->data.map.idx, e->data.map.mod, current_depth + e->data.map.depth,
-        current_row + e->data.map.row, e->data.map.col, e->buf, e->res_type);
+    int tid = 0;
+#if defined(_OPENMP)
+    tid = omp_get_thread_num();
+#endif
+    get_map_row(e->data.map.idx[tid], e->data.map.mod,
+                current_depth + e->data.map.depth,
+                current_row[tid] + e->data.map.row, e->data.map.col,
+                e->buf[tid], e->res_type);
 }
 
 static void evaluate_function(expression *e)
 {
     int i;
     int res;
+    int tid = 0;
+#if defined(_OPENMP)
+    tid = omp_get_thread_num();
+#endif
 
     if (e->data.func.argc > 1 && e->data.func.func != f_eval) {
         for (i = 1; i <= e->data.func.argc; i++)
@@ -216,8 +277,19 @@ static void evaluate_function(expression *e)
         for (i = 1; i <= e->data.func.argc; i++)
             evaluate(e->data.func.args[i]);
 
-    res = (*e->data.func.func)(e->data.func.argc, e->data.func.argt,
-                               e->data.func.argv);
+    /* copy the argv in the individual thread */
+    void **thread_argv = G_malloc((e->data.func.argc + 1) * sizeof(void *));
+    for (i = 0; i < e->data.func.argc + 1; i++)
+        thread_argv[i] = e->data.func.argv[i][tid];
+
+    res =
+        (*e->data.func.func)(e->data.func.argc, e->data.func.argt, thread_argv);
+
+    /* copy the results from thread_argv to e */
+    for (i = 0; i < e->data.func.argc + 1; i++)
+        e->data.func.argv[i][tid] = thread_argv[i];
+
+    G_free(thread_argv);
 
     switch (res) {
     case E_ARG_LO:
@@ -300,9 +372,12 @@ static void error_handler(void *p UNUSED)
 
 void execute(expr_list *ee)
 {
-    int verbose = isatty(2);
+    int verbose;
     expr_list *l;
-    int count, n;
+    expression **exp_arr;
+    int count, n, i;
+    int num_exprs = 0;
+    int threads = 1;
 
     exprs = ee;
     G_add_error_handler(error_handler, NULL);
@@ -323,13 +398,19 @@ void execute(expr_list *ee)
             G_fatal_error(_("output map <%s> exists. To overwrite, "
                             "use the --overwrite flag"),
                           var);
+        num_exprs++;
     }
 
+    /* Create a array of expreesion and stored it in heap */
+    exp_arr = G_malloc(num_exprs * sizeof(struct expression *));
+
     /* Parse each expression and extract all raster maps */
-    for (l = ee; l; l = l->next) {
+    l = ee;
+    for (i = 0; i < num_exprs; i++) {
         expression *e = l->exp;
-
         extract_maps(e);
+        exp_arr[i] = e;
+        l = l->next;
     }
 
     /* Set the region from the input maps */
@@ -341,8 +422,9 @@ void execute(expr_list *ee)
     setup_region();
 
     /* Parse each expression and initialize the maps, buffers and variables */
-    for (l = ee; l; l = l->next) {
-        expression *e = l->exp;
+
+    for (i = 0; i < num_exprs; i++) {
+        expression *e = exp_arr[i];
         const char *var;
         expression *val;
 
@@ -358,29 +440,50 @@ void execute(expr_list *ee)
 
     setup_maps();
 
+#if defined(_OPENMP)
+    threads = omp_get_max_threads();
+    /* Make sure the number of threads no more that the number of rows in
+     * rasters */
+    if ((threads > rows) && (threads > 1)) {
+        threads = rows;
+        omp_set_num_threads(threads);
+        G_verbose_message(
+            _("The number of rows is less than the number of threads. \
+            Set the number of threads to be the same as the rows = %d."),
+            threads);
+    }
+#endif
+    current_row = (int *)G_malloc(sizeof(int) * threads);
     count = rows * depths;
     n = 0;
 
-    G_init_workers();
-
+    verbose = isatty(2);
     for (current_depth = 0; current_depth < depths; current_depth++) {
-        for (current_row = 0; current_row < rows; current_row++) {
+#pragma omp parallel for default(shared) schedule(static, 1) private(i) ordered
+        for (int row = 0; row < rows; row++) {
             if (verbose)
                 G_percent(n, count, 2);
 
-            for (l = ee; l; l = l->next) {
-                expression *e = l->exp;
+            int tid = 0;
+#if defined(_OPENMP)
+            tid = omp_get_thread_num();
+#endif
+            /* calculate through expressions row by row */
+            current_row[tid] = row;
+            for (i = 0; i < num_exprs; i++) {
+                expression *e = exp_arr[i];
                 int fd;
-
                 evaluate(e);
-
-                if (e->type != expr_type_binding)
-                    continue;
-
-                fd = e->data.bind.fd;
-                put_map_row(fd, e->buf, e->res_type);
+#pragma omp ordered
+                {
+                    /* write out values to a file row by row */
+                    if (e->type == expr_type_binding) {
+                        fd = e->data.bind.fd;
+                        put_map_row(fd, e->buf[tid], e->res_type);
+                    }
+                }
             }
-
+#pragma omp atomic update
             n++;
         }
     }
@@ -410,17 +513,33 @@ void execute(expr_list *ee)
 
         if (val->type == expr_type_map) {
             if (val->data.map.mod == 'M') {
-                copy_cats(var, val->data.map.idx);
-                copy_colors(var, val->data.map.idx);
+                copy_cats(var, val->data.map.idx[0]);
+                copy_colors(var, val->data.map.idx[0]);
             }
 
-            copy_history(var, val->data.map.idx);
+            copy_history(var, val->data.map.idx[0]);
         }
         else
             create_history(var, val);
     }
 
     G_unset_error_routine();
+
+    /* Free the memory and make it unreachable */
+    G_free(current_row);
+    for (i = 0; i < num_exprs; i++) {
+        expression *e = exp_arr[i];
+        free_buf(e);
+        if (e->type == expr_type_function)
+            free_argv(e);
+        if (e->type == expr_type_map && e->data.map.idx) {
+            G_free(e->data.map.idx);
+            e->data.map.idx = NULL;
+        }
+    }
+    G_free(exp_arr);
+    current_row = NULL;
+    exp_arr = NULL;
 }
 
 void describe_maps(FILE *fp, expr_list *ee)
diff --git a/raster/r.mapcalc/expression.h b/raster/r.mapcalc/expression.h
index 9f049cc7a6a..f55fffd3650 100644
--- a/raster/r.mapcalc/expression.h
+++ b/raster/r.mapcalc/expression.h
@@ -27,7 +27,7 @@ typedef struct expr_data_map {
     const char *name;
     int mod;
     int row, col, depth;
-    int idx;
+    int *idx; /* array to store fds for multi-threads*/
 } expr_data_map;
 
 typedef struct expr_data_func {
@@ -35,10 +35,10 @@ typedef struct expr_data_func {
     const char *oper;
     int prec;
     func_t *func;
-    int argc;
-    struct expression **args;
-    int *argt;
-    void **argv;
+    int argc;                 /* number of args in the whole expression */
+    struct expression **args; /* array of expressions */
+    int *argt;                /* type of expressions */
+    void ***argv;             /* values in e->buf for each expression */
 } expr_data_func;
 
 typedef struct expr_data_bind {
@@ -50,7 +50,7 @@ typedef struct expr_data_bind {
 typedef struct expression {
     int type;
     int res_type;
-    void *buf;
+    void **buf;
     union {
         expr_data_const con;
         expr_data_var var;
diff --git a/raster/r.mapcalc/globals.h b/raster/r.mapcalc/globals.h
index 0d49c166198..351ebbabe36 100644
--- a/raster/r.mapcalc/globals.h
+++ b/raster/r.mapcalc/globals.h
@@ -6,7 +6,8 @@ extern long seed_value;
 extern long seeded;
 extern int region_approach;
 
-extern int current_depth, current_row;
+extern int current_depth;
+extern int *current_row;
 extern int depths, rows, columns;
 
 #endif /* __GLOBALS_H_ */
diff --git a/raster/r.mapcalc/main.c b/raster/r.mapcalc/main.c
index e18214604ae..faa2fed12c1 100644
--- a/raster/r.mapcalc/main.c
+++ b/raster/r.mapcalc/main.c
@@ -11,6 +11,9 @@
  *               for details.
  *
  *****************************************************************************/
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
 
 #include <unistd.h>
 #include <stdio.h>
@@ -59,10 +62,11 @@ static expr_list *parse_file(const char *filename)
 int main(int argc, char **argv)
 {
     struct GModule *module;
-    struct Option *expr, *file, *seed, *region;
+    struct Option *expr, *file, *seed, *region, *nprocs;
     struct Flag *random, *describe;
     int all_ok;
     char *desc;
+    int threads = 1;
 
     G_gisinit(argv[0]);
 
@@ -115,9 +119,10 @@ int main(int argc, char **argv)
     describe->key = 'l';
     describe->description = _("List input and output maps");
 
-    if (argc == 1) {
-        char **p = G_malloc(3 * sizeof(char *));
+    nprocs = G_define_standard_option(G_OPT_M_NPROCS);
 
+    char **p = G_malloc(3 * sizeof(char *));
+    if (argc == 1) {
         p[0] = argv[0];
         p[1] = G_store("file=-");
         p[2] = NULL;
@@ -179,11 +184,45 @@ int main(int argc, char **argv)
     }
 
     pre_exec();
+
+    /* Determine the number of threads */
+    threads = atoi(nprocs->answer);
+
+    /* Check if the program name is r3.mapcalc */
+    /* Handle both Unix and Windows path separators */
+    const char *progname = strrchr(argv[0], '/');
+    if (!progname)
+        progname = strrchr(argv[0], '\\');
+    progname = progname ? progname + 1 : argv[0];
+
+    if ((strncmp(progname, "r3.mapcalc", 10) == 0) && (threads != 1)) {
+        threads = 1;
+        nprocs->answer = "1";
+        G_verbose_message(_("r3.mapcalc does not support parallel execution."));
+    }
+    else if ((seeded) && (threads != 1)) {
+        threads = 1;
+        nprocs->answer = "1";
+        G_verbose_message(
+            _("Parallel execution is not supported for random seed."));
+    }
+
+    /* Ensure the proper number of threads is assigned */
+    threads = G_set_omp_num_threads(nprocs);
+    if (threads > 1)
+        threads = Rast_disable_omp_on_mask(threads);
+    if (threads < 1)
+        G_fatal_error(_("<%d> is not valid number of nprocs."), threads);
+
+    /* Execute calculations */
     execute(result);
     post_exec();
 
     all_ok = 1;
 
+    G_free(p);
+    p = NULL;
+
     if (floating_point_exception_occurred) {
         G_warning(_("Floating point error(s) occurred in the calculation"));
         all_ok = 0;
diff --git a/raster/r.mapcalc/map.c b/raster/r.mapcalc/map.c
index 617b8f358a3..bb6475a1f6e 100644
--- a/raster/r.mapcalc/map.c
+++ b/raster/r.mapcalc/map.c
@@ -1,3 +1,7 @@
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
+
 #include <grass/config.h>
 
 #include <stdlib.h>
@@ -59,7 +63,7 @@ struct map {
     struct Categories cats;
     struct Colors colors;
     BTREE btree;
-    struct row_cache cache;
+    struct row_cache *caches;
 #ifdef HAVE_PTHREAD_H
     pthread_mutex_t mutex;
 #endif
@@ -148,10 +152,7 @@ static void cache_release(struct row_cache *cache)
 static void *cache_get_raw(struct row_cache *cache, int row, int data_type)
 {
     struct sub_cache *sub;
-    void **tmp;
-    char *vtmp;
-    int i, j;
-    int newrow;
+    int i;
 
     if (!cache->sub[data_type])
         cache_sub_init(cache, data_type);
@@ -175,31 +176,11 @@ static void *cache_get_raw(struct row_cache *cache, int row, int data_type)
         return sub->buf[0];
     }
 
-    tmp = G_alloca(cache->nrows * sizeof(void *));
-    memcpy(tmp, sub->buf, cache->nrows * sizeof(void *));
-    vtmp = G_alloca(cache->nrows);
-    memcpy(vtmp, sub->valid, cache->nrows);
-
-    i = (i < 0) ? 0 : cache->nrows - 1;
-    newrow = row - i;
-
-    for (j = 0; j < cache->nrows; j++) {
-        int r = newrow + j;
-        int k = r - sub->row;
-        int l = (k + cache->nrows) % cache->nrows;
-
-        sub->buf[j] = tmp[l];
-        sub->valid[j] = k >= 0 && k < cache->nrows && vtmp[l];
+    else {
+        i = (i < 0) ? 0 : cache->nrows - 1;
+        read_row(cache->fd, sub->buf[i], row, data_type);
+        return sub->buf[i];
     }
-
-    sub->row = newrow;
-    G_freea(tmp);
-    G_freea(vtmp);
-
-    read_row(cache->fd, sub->buf[i], row, data_type);
-    sub->valid[i] = 1;
-
-    return sub->buf[i];
 }
 
 static void cache_get(struct row_cache *cache, void *buf, int row, int res_type)
@@ -382,13 +363,19 @@ static void translate_from_cats(struct map *m, CELL *cell, DCELL *xcell,
 static void setup_map(struct map *m)
 {
     int nrows = m->max_row - m->min_row + 1;
-
+    int threads = 1;
 #ifdef HAVE_PTHREAD_H
     pthread_mutex_init(&m->mutex, NULL);
 #endif
 
+#ifdef _OPENMP
+    threads = omp_get_max_threads();
+#endif
+    m->caches =
+        (struct row_cache *)G_malloc(threads * sizeof(struct row_cache));
     if (nrows > 1 && nrows <= max_rows_in_memory) {
-        cache_setup(&m->cache, m->fd, nrows);
+        for (int i = 0; i < threads; i++)
+            cache_setup(&m->caches[i], m->fd, nrows);
         m->use_rowio = 1;
     }
     else
@@ -425,17 +412,30 @@ static void read_map(struct map *m, void *buf, int res_type, int row, int col)
         return;
     }
 
+    int tid = 0;
+#ifdef _OPENMP
+    tid = omp_get_thread_num();
+#endif
+
     if (m->use_rowio)
-        cache_get(&m->cache, buf, row, res_type);
+        cache_get(&m->caches[tid], buf, row, res_type);
     else
         read_row(m->fd, buf, row, res_type);
 
+#ifdef _OPENMP
+    tid = omp_get_thread_num();
+#endif
     if (col)
         column_shift(buf, res_type, col);
 }
 
 static void close_map(struct map *m)
 {
+    int threads = 1;
+#ifdef _OPENMP
+    threads = omp_get_max_threads();
+#endif
+
     if (m->fd < 0)
         return;
 
@@ -457,7 +457,10 @@ static void close_map(struct map *m)
     }
 
     if (m->use_rowio) {
-        cache_release(&m->cache);
+        for (int i = 0; i < threads; i++)
+            cache_release(&m->caches[i]);
+        if (threads > 1)
+            G_free(m->caches);
         m->use_rowio = 0;
     }
 }
@@ -494,7 +497,6 @@ int map_type(const char *name, int mod)
 
 int open_map(const char *name, int mod, int row, int col)
 {
-    int i;
     const char *mapset;
     int use_cats = 0;
     int use_colors = 0;
@@ -532,26 +534,6 @@ int open_map(const char *name, int mod, int row, int col)
         break;
     }
 
-    for (i = 0; i < num_maps; i++) {
-        m = &maps[i];
-
-        if (strcmp(m->name, name) != 0 || strcmp(m->mapset, mapset) != 0)
-            continue;
-
-        if (row < m->min_row)
-            m->min_row = row;
-        if (row > m->max_row)
-            m->max_row = row;
-
-        if (use_cats && !m->have_cats)
-            init_cats(m);
-
-        if (use_colors && !m->have_colors)
-            init_colors(m);
-
-        return i;
-    }
-
     if (num_maps >= max_maps) {
         max_maps += 10;
         maps = G_realloc(maps, max_maps * sizeof(struct map));
diff --git a/raster/r.mapcalc/map3.c b/raster/r.mapcalc/map3.c
index 7bde2a58cce..9a4dbd20b48 100644
--- a/raster/r.mapcalc/map3.c
+++ b/raster/r.mapcalc/map3.c
@@ -1,3 +1,7 @@
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
+
 #include <stdlib.h>
 #include <limits.h>
 #include <string.h>
@@ -619,8 +623,12 @@ int open_output_map(const char *name, int res_type)
 void put_map_row(int fd, void *buf, int res_type)
 {
     void *handle = omaps[fd];
+    int tid = 0;
+#if defined(_OPENMP)
+    tid = omp_get_thread_num();
+#endif
 
-    write_row(handle, buf, res_type, current_depth, current_row);
+    write_row(handle, buf, res_type, current_depth, current_row[tid]);
 }
 
 void close_output_map(int fd)
diff --git a/raster/r.mapcalc/mapcalc.h b/raster/r.mapcalc/mapcalc.h
index 75d6b3da020..a2c6947c641 100644
--- a/raster/r.mapcalc/mapcalc.h
+++ b/raster/r.mapcalc/mapcalc.h
@@ -2,6 +2,9 @@
 #define _MAPCALC_H_
 
 /****************************************************************************/
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
 
 #include <stdio.h>
 
diff --git a/raster/r.mapcalc/r.mapcalc.md b/raster/r.mapcalc/r.mapcalc.md
index fd808d8c384..e2766668f56 100644
--- a/raster/r.mapcalc/r.mapcalc.md
+++ b/raster/r.mapcalc/r.mapcalc.md
@@ -872,6 +872,26 @@ X (map) values supplied and y (newmap) values returned:
 100, 50
 ```
 
+### Performance
+
+r.mapcalc is parallelized using OpenMP. The number of threads can be controlled
+with the **nprocs** parameter. Note that more complex expressions can benefit
+from more threads. By default (**nprocs=0**), r.mapcalc uses all available threads.
+Use the **--verbose** flag to display the number of threads in use.
+If you observe reduced performance when using many threads, try lowering ther number.
+
+Note: r.mapcalc may disable parallelization in certain cases, even when requested:
+
+- When a mask is active, because the current parallel implementation
+ does not support it.
+- When the rand() function is used, to ensure reproducible results.
+
+![Benchmark of r.mapcalc](r_mapcalc_benchmark_time.png)  
+*Figure: Benchmark shows execution time for different number of cells
+and different complexity of expressions.
+See benchmark script in the source code.
+(Intel Core i9-10940X CPU @ 3.30GHz x 28)*
+
 ## KNOWN ISSUES
 
 The *result* variable on the left hand side of the equation should not
diff --git a/raster/r.mapcalc/r_mapcalc_benchmark_time.png b/raster/r.mapcalc/r_mapcalc_benchmark_time.png
new file mode 100644
index 00000000000..c499d491fda
Binary files /dev/null and b/raster/r.mapcalc/r_mapcalc_benchmark_time.png differ
diff --git a/raster/r.mapcalc/testsuite/test_nmedian_bug_3296.py b/raster/r.mapcalc/testsuite/test_nmedian_bug_3296.py
index 3ab7d7c255f..c5a8e8f4690 100644
--- a/raster/r.mapcalc/testsuite/test_nmedian_bug_3296.py
+++ b/raster/r.mapcalc/testsuite/test_nmedian_bug_3296.py
@@ -125,6 +125,41 @@ def test_dcell(self):
             actual=self.output, reference=self.output_ref, precision=0
         )
 
+    def test_cell_nprocs1(self):
+        expression = "{o}=nmedian(({i}[0,-1] - {i})^2,({i}[0,1] - {i})^2)".format(
+            o=self.output, i=self.input
+        )
+        self.assertModule("r.mapcalc", expression=expression, nprocs=1, overwrite=True)
+        self.assertRasterExists(self.output)
+        self.to_remove.append(self.output)
+        self.assertRastersNoDifference(
+            actual=self.output, reference=self.output_cell, precision=0
+        )
+
+    def test_fcell_nprocs1(self):
+        expression = (
+            "{o}=nmedian(float(({i}[0,-1] - {i})^2), float(({i}[0,1] - {i})^2))".format(
+                o=self.output, i=self.input
+            )
+        )
+        self.assertModule("r.mapcalc", expression=expression, nprocs=1, overwrite=True)
+        self.assertRasterExists(self.output)
+        self.to_remove.append(self.output)
+        self.assertRastersNoDifference(
+            actual=self.output, reference=self.output_ref, precision=0
+        )
+
+    def test_dcell_nprocs1(self):
+        expression = "{o}=nmedian(double(({i}[0,-1] - {i})^2), double(({i}[0,1] - {i})^2))".format(
+            o=self.output, i=self.input
+        )
+        self.assertModule("r.mapcalc", expression=expression, nprocs=1, overwrite=True)
+        self.assertRasterExists(self.output)
+        self.to_remove.append(self.output)
+        self.assertRastersNoDifference(
+            actual=self.output, reference=self.output_ref, precision=0
+        )
+
 
 if __name__ == "__main__":
     test()
diff --git a/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py b/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py
new file mode 100644
index 00000000000..c635785daba
--- /dev/null
+++ b/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py
@@ -0,0 +1,367 @@
+from grass.gunittest.case import TestCase
+from grass.gunittest.main import test
+from grass.gunittest.gmodules import SimpleModule
+
+cell_seed_500 = """\
+north: 20
+south: 10
+east: 25
+west: 15
+rows: 10
+cols: 10
+121 12 183 55 37 96 138 117 182 40
+157 70 115 1 149 125 42 193 108 24
+83 66 82 84 186 182 179 122 67 113
+151 93 144 173 128 196 61 125 64 193
+180 175 14 41 44 27 165 27 90 60
+97 57 12 104 98 13 87 24 83 107
+174 133 146 114 115 60 78 154 49 130
+55 138 144 25 32 58 47 137 139 32
+143 193 155 190 131 124 87 81 160 154
+56 45 48 66 9 182 69 12 154 19
+"""
+
+dcell_seed_600 = """\
+north: 20
+south: 10
+east: 25
+west: 15
+rows: 10
+cols: 10
+130.790433856418332 101.3319248101491041 33.5781271447787759 37.4064724824657944 98.2794723130458152 73.9118866262841863 185.9530433718733775 74.5210037729812882 166.1178416001017695 90.9915650902159996
+109.2478664232956334 25.6499350759712215 150.9024447059825036 125.7544119036241312 66.7235333366722614 167.9375729129454271 123.1009291055983965 12.0922254083554606 59.389026967287819 113.2843489528100349
+40.0044184023145277 135.8273774212801186 71.6737798852435049 191.6223505280372876 4.1546013811569615 143.3082794522489962 177.043829294835632 115.0300571162354402 141.8985452774071234 127.8949967123638061
+93.2842559637482793 9.7471880052423856 118.1216452002055632 158.1474162140586088 67.2957262519499437 3.6524546350146849 147.0965842667525862 37.060628529871579 47.3408278816968959 66.2219633495724054
+175.5638637866295539 67.1399023507611901 162.2058392782793703 198.1586789345953719 36.474049475167746 49.2589028048889617 112.1169663235969836 22.0227597984432535 95.9169228571662131 86.7470895014531322
+93.5401613888204935 193.7821104138942587 193.8286564351004699 3.2623643889134684 94.6955247357847725 25.7099391122614307 155.592251526775442 25.3392337002970294 48.3979699868663005 99.6836079482556272
+104.16296861365457 190.7865884377180805 6.2841805474238619 49.3731395705159528 100.1903962703459285 116.927654961282343 19.8626348109264264 40.9693022766258466 81.6500759554420057 169.2220572316770131
+118.8112518721558217 55.8955021401724039 112.9150308215961331 62.6399760484719081 85.400498505854145 191.0144187084912062 124.2128169358724534 167.9341741649760706 170.6149695243870781 158.3034517206661462
+130.0453795775294736 64.1996403829061535 62.9317494959142465 175.1909990236256931 122.9624852869890361 79.9546265736285733 9.6594013716963367 114.0611338072915544 11.9371167643030809 186.9121199748369122
+3.2891990250261536 30.9245408751958379 46.4021422454598991 104.2378950097200203 47.424093232347019 73.4801303522840499 22.4778583078695213 132.870185207462697 48.1666164169167388 100.5504714442693057
+"""
+
+fcell_seed_700 = """\
+north: 20
+south: 10
+east: 25
+west: 15
+rows: 10
+cols: 10
+146.756378 192.682159 2.644822 147.270462 62.178818 192.668198 94.320778 107.710426 98.319664 114.444504
+12.995321 18.026272 151.590958 5.249451 197.266708 103.663635 115.424088 28.01062 78.555168 62.912098
+164.053619 154.652039 98.536011 44.601639 85.322289 168.383957 44.93845 128.62262 89.910591 107.242188
+111.182487 63.080284 177.791473 47.439354 42.451859 72.396568 170.597778 170.622742 141.88858 105.126854
+120.76828 148.581085 42.124866 56.432236 164.652176 98.094009 60.741329 66.286987 187.847427 160.120056
+50.530689 179.090652 138.114014 138.629211 193.147903 172.861481 133.72728 108.720459 103.508438 28.81559
+39.653179 101.948265 35.744762 25.570076 78.767021 154.600616 144.907684 82.370148 116.378654 18.218494
+35.587288 66.534409 65.744408 186.476959 137.081116 151.379272 48.261463 8.323328 130.432739 53.346546
+152.67189 15.512391 146.049072 185.276245 34.417141 127.522453 124.54998 52.08218 167.141342 87.771118
+69.0522 43.57811 63.15279 68.677063 74.202805 97.429077 167.123199 19.892767 120.593437 190.960815
+"""
+
+THREADS = 4
+
+
+class TestRandFunction(TestCase):
+    # TODO: replace by unified handing of maps
+    to_remove = []
+
+    @classmethod
+    def setUpClass(cls):
+        cls.use_temp_region()
+        cls.runModule("g.region", n=20, s=10, e=25, w=15, res=1)
+
+    @classmethod
+    def tearDownClass(cls):
+        cls.del_temp_region()
+        if cls.to_remove:
+            cls.runModule(
+                "g.remove", flags="f", type="raster", name=",".join(cls.to_remove)
+            )
+
+    def rinfo_contains_number(self, raster, number):
+        """Test that r.info standard output for raster contains a given number
+
+        To be used in test methods for testing presence of a given number.
+        """
+        rinfo = SimpleModule("r.info", map=raster)
+        self.runModule(rinfo)
+        self.assertIn(str(number), rinfo.outputs.stdout)
+
+    def test_seed_not_required(self):
+        """Test that seed is not required when rand() is not used"""
+        self.assertModule("r.mapcalc", expression="nonrand_cell = 200", nprocs=THREADS)
+        self.to_remove.append("nonrand_cell")
+
+    def test_seed_required(self):
+        """Test that seed is required when rand() is used
+
+        This test can, and probably should, generate an error message.
+        """
+        self.assertModuleFail(
+            "r.mapcalc", expression="rand_x = rand(1, 200)", nprocs=THREADS
+        )
+        # TODO: assert map not exists but it would be handy here
+        # TODO: test that error message was generated
+
+    def test_seed_cell(self):
+        """Test given seed with CELL against reference map"""
+        seed = 500
+        self.runModule(
+            "r.in.ascii", input="-", stdin=cell_seed_500, output="rand_cell_ref"
+        )
+        self.to_remove.append("rand_cell_ref")
+        self.assertModule(
+            "r.mapcalc",
+            seed=seed,
+            expression="rand_cell = rand(1, 200)",
+            nprocs=THREADS,
+        )
+        self.to_remove.append("rand_cell")
+        # this assert is using r.mapcalc but we are testing different
+        # functionality than used by assert
+        self.assertRastersNoDifference(
+            actual="rand_cell", reference="rand_cell_ref", precision=0
+        )
+        self.rinfo_contains_number("rand_cell", seed)
+
+    def test_seed_dcell(self):
+        """Test given seed with DCELL against reference map"""
+        seed = 600
+        self.runModule(
+            "r.in.ascii", input="-", stdin=dcell_seed_600, output="rand_dcell_ref"
+        )
+        self.to_remove.append("rand_dcell_ref")
+        self.assertModule(
+            "r.mapcalc",
+            seed=seed,
+            expression="rand_dcell = rand(1.0, 200.0)",
+            nprocs=THREADS,
+        )
+        self.to_remove.append("rand_dcell")
+        # this assert is using r.mapcalc but we are testing different
+        # functionality than used by assert
+        self.assertRastersNoDifference(
+            actual="rand_dcell", reference="rand_dcell_ref", precision=0.00000000000001
+        )
+        self.rinfo_contains_number("rand_dcell", seed)
+
+    def test_seed_fcell(self):
+        """Test given seed with FCELL against reference map"""
+        seed = 700
+        self.runModule(
+            "r.in.ascii", input="-", stdin=fcell_seed_700, output="rand_fcell_ref"
+        )
+        self.to_remove.append("rand_fcell_ref")
+        self.assertModule(
+            "r.mapcalc",
+            seed=seed,
+            expression="rand_fcell = rand(float(1), 200)",
+            nprocs=THREADS,
+        )
+        self.to_remove.append("rand_fcell")
+        # this assert is using r.mapcalc but we are testing different
+        # functionality than used by assert
+        self.assertRastersNoDifference(
+            actual="rand_fcell", reference="rand_fcell_ref", precision=0.000001
+        )
+        self.rinfo_contains_number("rand_fcell", seed)
+
+    def test_auto_seed(self):
+        """Test that two runs with -s does not give same maps"""
+        self.assertModule(
+            "r.mapcalc",
+            flags="s",
+            expression="rand_auto_1 = rand(1., 2)",
+            nprocs=THREADS,
+        )
+        self.to_remove.append("rand_auto_1")
+        self.assertModule(
+            "r.mapcalc",
+            flags="s",
+            expression="rand_auto_2 = rand(1., 2)",
+            nprocs=THREADS,
+        )
+        self.to_remove.append("rand_auto_2")
+        self.assertRastersDifference(
+            "rand_auto_1",
+            "rand_auto_2",
+            statistics={"min": -1, "max": 1, "mean": 0},
+            precision=0.5,
+        )  # low precision, we have few cells
+
+
+# TODO: add more expressions
+# TODO: add tests with prepared data
+
+
+class TestBasicOperations(TestCase):
+    # TODO: replace by unified handing of maps
+    to_remove = []
+
+    @classmethod
+    def setUpClass(cls):
+        cls.use_temp_region()
+        cls.runModule("g.region", n=20, s=10, e=25, w=15, res=1)
+
+    @classmethod
+    def tearDownClass(cls):
+        cls.del_temp_region()
+        if cls.to_remove:
+            cls.runModule(
+                "g.remove",
+                flags="f",
+                type="raster",
+                name=",".join(cls.to_remove),
+                verbose=True,
+            )
+
+    def test_difference_of_the_same_map_double(self):
+        """Test zero difference of map with itself"""
+        self.runModule("r.mapcalc", flags="s", expression="a = rand(1.0, 200)")
+        self.to_remove.append("a")
+        self.assertModule("r.mapcalc", expression="diff_a_a = a - a", nprocs=THREADS)
+        self.to_remove.append("diff_a_a")
+        self.assertRasterMinMax("diff_a_a", refmin=0, refmax=0)
+
+    def test_difference_of_the_same_map_float(self):
+        """Test zero difference of map with itself"""
+        self.runModule("r.mapcalc", flags="s", expression="af = rand(float(1), 200)")
+        self.to_remove.append("af")
+        self.assertModule(
+            "r.mapcalc", expression="diff_af_af = af - af", nprocs=THREADS
+        )
+        self.to_remove.append("diff_af_af")
+        self.assertRasterMinMax("diff_af_af", refmin=0, refmax=0)
+
+    def test_difference_of_the_same_map_int(self):
+        """Test zero difference of map with itself"""
+        self.runModule("r.mapcalc", flags="s", expression="ai = rand(1, 200)")
+        self.to_remove.append("ai")
+        self.assertModule(
+            "r.mapcalc", expression="diff_ai_ai = ai - ai", nprocs=THREADS
+        )
+        self.to_remove.append("diff_ai_ai")
+        self.assertRasterMinMax("diff_ai_ai", refmin=0, refmax=0)
+
+    def test_difference_of_the_same_expression(self):
+        """Test zero difference of two same expressions"""
+        self.assertModule(
+            "r.mapcalc",
+            expression="diff_e_e = 3 * x() * y() - 3 * x() * y()",
+            nprocs=THREADS,
+        )
+        self.to_remove.append("diff_e_e")
+        self.assertRasterMinMax("diff_e_e", refmin=0, refmax=0)
+
+    def test_nrows_ncols_sum(self):
+        """Test if sum of nrows and ncols matches one
+        expected from current region settings"""
+        self.assertModule(
+            "r.mapcalc",
+            expression="nrows_ncols_sum = nrows() + ncols()",
+            nprocs=THREADS,
+        )
+        self.to_remove.append("nrows_ncols_sum")
+        self.assertRasterMinMax("nrows_ncols_sum", refmin=20, refmax=20)
+
+
+class TestRegionOperations(TestCase):
+    # TODO: replace by unified handing of maps
+    to_remove = []
+
+    @classmethod
+    def setUpClass(cls):
+        cls.use_temp_region()
+        cls.runModule("g.region", n=30, s=15, e=30, w=15, res=5)
+        cls.runModule(
+            "r.mapcalc", expression="test_region_1 = 1", seed=1, nprocs=THREADS
+        )
+        cls.runModule("g.region", n=25, s=10, e=25, w=10, res=5)
+        cls.runModule(
+            "r.mapcalc", expression="test_region_2 = 2", seed=1, nprocs=THREADS
+        )
+        cls.runModule("g.region", n=20, s=5, e=20, w=5, res=1)
+        cls.runModule(
+            "r.mapcalc", expression="test_region_3 = 3", seed=1, nprocs=THREADS
+        )
+
+        cls.to_remove.append("test_region_1")
+        cls.to_remove.append("test_region_2")
+        cls.to_remove.append("test_region_3")
+
+    @classmethod
+    def tearDownClass(cls):
+        cls.del_temp_region()
+        if cls.to_remove:
+            cls.runModule(
+                "g.remove",
+                flags="f",
+                type="raster",
+                name=",".join(cls.to_remove),
+                verbose=True,
+            )
+
+    def test_union(self):
+        """Test the union region option"""
+        self.assertModule(
+            "r.mapcalc",
+            region="union",
+            seed=1,
+            expression="test_region_4 = test_region_1 + test_region_2 + test_region_3",
+            nprocs=THREADS,
+        )
+        self.to_remove.append("test_region_4")
+
+        self.assertModuleKeyValue(
+            "r.info",
+            map="test_region_4",
+            flags="gr",
+            reference={
+                "min": 6,
+                "max": 6,
+                "cells": 625,
+                "north": 30,
+                "south": 5,
+                "west": 5,
+                "east": 30,
+                "nsres": 1,
+                "ewres": 1,
+            },
+            precision=0.01,
+            sep="=",
+        )
+
+    def test_intersect(self):
+        """Test the intersect region option"""
+        self.assertModule(
+            "r.mapcalc",
+            region="intersect",
+            seed=1,
+            expression="test_region_5 = test_region_1 + test_region_2 + test_region_3",
+            nprocs=THREADS,
+        )
+        self.to_remove.append("test_region_5")
+
+        self.assertModuleKeyValue(
+            "r.info",
+            map="test_region_5",
+            flags="gr",
+            reference={
+                "min": 6,
+                "max": 6,
+                "cells": 25,
+                "north": 20,
+                "south": 15,
+                "west": 15,
+                "east": 20,
+                "nsres": 1,
+                "ewres": 1,
+            },
+            precision=0.01,
+            sep="=",
+        )
+
+
+if __name__ == "__main__":
+    test()
diff --git a/raster/r.mapcalc/xarea.c b/raster/r.mapcalc/xarea.c
index 74bd1392c3c..6bbe605ba6a 100644
--- a/raster/r.mapcalc/xarea.c
+++ b/raster/r.mapcalc/xarea.c
@@ -1,3 +1,7 @@
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
+
 #include <grass/gis.h>
 #include <grass/raster.h>
 #include "globals.h"
@@ -10,6 +14,11 @@ area() area of a cell in square meters
 
 int f_area(int argc, const int *argt, void **args)
 {
+    int tid = 0;
+#if defined(_OPENMP)
+    tid = omp_get_thread_num();
+#endif
+
     DCELL *res = args[0];
     int i;
     static int row = -1;
@@ -21,11 +30,11 @@ int f_area(int argc, const int *argt, void **args)
     if (argt[0] != DCELL_TYPE)
         return E_RES_TYPE;
 
-    if (row != current_row) {
+    if (row != current_row[tid]) {
         if (row == -1)
             G_begin_cell_area_calculations();
 
-        row = current_row;
+        row = current_row[tid];
         cell_area = G_area_of_cell_at_row(row);
     }
 
diff --git a/raster/r.mapcalc/xcoor.c b/raster/r.mapcalc/xcoor.c
index 5f783d7e90f..1a209a25281 100644
--- a/raster/r.mapcalc/xcoor.c
+++ b/raster/r.mapcalc/xcoor.c
@@ -1,3 +1,7 @@
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
+
 #include <grass/gis.h>
 #include <grass/raster.h>
 #include "globals.h"
@@ -35,6 +39,11 @@ int f_x(int argc, const int *argt, void **args)
 
 int f_y(int argc, const int *argt, void **args)
 {
+    int tid = 0;
+#if defined(_OPENMP)
+    tid = omp_get_thread_num();
+#endif
+
     DCELL *res = args[0];
     DCELL y;
     int i;
@@ -45,7 +54,7 @@ int f_y(int argc, const int *argt, void **args)
     if (argt[0] != DCELL_TYPE)
         return E_RES_TYPE;
 
-    y = Rast_row_to_northing(current_row + 0.5, &current_region2);
+    y = Rast_row_to_northing(current_row[tid] + 0.5, &current_region2);
 
     for (i = 0; i < columns; i++)
         res[i] = y;
diff --git a/raster/r.mapcalc/xcoor3.c b/raster/r.mapcalc/xcoor3.c
index 8c8906b890e..21fa5a381b0 100644
--- a/raster/r.mapcalc/xcoor3.c
+++ b/raster/r.mapcalc/xcoor3.c
@@ -1,3 +1,7 @@
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
+
 #include <grass/gis.h>
 #include <grass/raster3d.h>
 #include "globals.h"
@@ -36,6 +40,11 @@ int f_x(int argc, const int *argt, void **args)
 
 int f_y(int argc, const int *argt, void **args)
 {
+    int tid = 0;
+#if defined(_OPENMP)
+    tid = omp_get_thread_num();
+#endif
+
     RASTER3D_Region *window = &current_region3;
     DCELL *res = args[0];
     DCELL y;
@@ -47,7 +56,7 @@ int f_y(int argc, const int *argt, void **args)
     if (argt[0] != DCELL_TYPE)
         return E_RES_TYPE;
 
-    y = window->north - (current_row + 0.5) * window->ns_res;
+    y = window->north - (current_row[tid] + 0.5) * window->ns_res;
 
     for (i = 0; i < columns; i++)
         res[i] = y;
diff --git a/raster/r.mapcalc/xrowcol.c b/raster/r.mapcalc/xrowcol.c
index 5f95d63e618..4f332589e79 100644
--- a/raster/r.mapcalc/xrowcol.c
+++ b/raster/r.mapcalc/xrowcol.c
@@ -1,3 +1,7 @@
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
+
 #include <grass/gis.h>
 #include <grass/raster.h>
 #include "globals.h"
@@ -32,8 +36,13 @@ int f_col(int argc, const int *argt, void **args)
 
 int f_row(int argc, const int *argt, void **args)
 {
+    int tid = 0;
+#if defined(_OPENMP)
+    tid = omp_get_thread_num();
+#endif
+
     CELL *res = args[0];
-    int row = current_row + 1;
+    int row = current_row[tid] + 1;
     int i;
 
     if (argc > 0)