r.sun: port r.sun.mp by Jaro Hofierka et al. from addons to core and add tests

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@71611 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.sun/Makefile

@@ -2,9 +2,10 @@ MODULE_TOPDIR = ../..
 
 PGM = r.sun
 
-LIBES = $(GPROJLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB) $(PROJLIB)
-DEPENDENCIES = $(GPROJDEP) $(RASTERDEP) $(GISDEP)
+LIBES = $(GPROJLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB) $(GMATHLIB) $(PROJLIB) $(OMPLIB)
+DEPENDENCIES = $(GPROJDEP) $(RASTERDEP) $(GISDEP) $(GMATHDEP)
 EXTRA_INC = $(PROJINC) $(OCLINCPATH)
+EXTRA_CFLAGS = $(OMPCFLAGS)
 EXTRA_LIBS = $(OCLLIB)
 #needed?  LIBES += $(OCLLIBPATH)  or  EXTRA_LDFLAGS = $(OCLLIBPATH)
 

raster/r.sun/local_proto.h

@@ -3,8 +3,31 @@
 
 void where_is_point(double *length, struct SunGeometryVarDay *sunVarGeom,
 		    struct GridGeometry *gridGeom);
+void where_is_point_parallel(double *length, double *sunVarGeom_zp,
+                             double *gridGeom_xx0,
+                             double *gridGeom_yy0,
+                             double *gridGeom_xg0,
+                             double *gridGeom_yg0,
+                             double *gridGeom_stepx,
+                             double *gridGeom_stepy);
 int searching(double *length, struct SunGeometryVarDay *sunVarGeom,
 	      struct GridGeometry *gridGeom);
+int searching_parallel(double *length,
+                       double *sunVarGeom_z_orig,
+                       double *sunVarGeom_zmax,
+                       double *sunVarGeom_zp,
+                       double *sunVarGeom_tanSolarAltitude,
+                       double *sunVarGeom_stepsinangle,
+                       double *sunVarGeom_stepcosangle,
+                       double *gridGeom_xx0,
+                       double *gridGeom_yy0,
+                       double *gridGeom_xg0,
+                       double *gridGeom_yg0,
+                       double *gridGeom_stepx,
+                       double *gridGeom_stepy,
+                       double *gridGeom_deltx,
+                       double *gridGeom_delty);
+
 
 int useCivilTime();
 void setUseCivilTime(int val);
@@ -29,11 +52,31 @@ double com_sol_const(int no_of_day);
 double brad(double, double *bh, struct SunGeometryVarDay *sunVarGeom,
 	    struct SunGeometryVarSlope *sunSlopeGeom,
 	    struct SolarRadVar *sunRadVar);
+double brad_parallel(double sh, double *bh,
+                     double *sunVarGeom_z_orig,
+                     double *sunVarGeom_solarAltitude,
+                     double *sunVarGeom_sinSolarAltitude,
+                     double *sunSlopeGeom_slope,
+                     double *sunSlopeGeom_aspect,
+                     double *sunRadVar_cbh,
+                     double *sunRadVar_linke,
+                     double *sunRadVar_G_norm_extra);
+
 double drad(double, double bh, double *rr,
 	    struct SunGeometryVarDay *sunVarGeom,
 	    struct SunGeometryVarSlope *sunSlopeGeom,
 	    struct SolarRadVar *sunRadVar);
-
+double drad_parallel(double sh, double bh, double *rr,
+                     int *sunVarGeom_isShadow,
+                     double *sunVarGeom_solarAltitude,
+                     double *sunVarGeom_sinSolarAltitude,
+                     double *sunVarGeom_solarAzimuth,
+                     double *sunSlopeGeom_aspect,
+                     double *sunSlopeGeom_slope,
+                     double *sunRadVar_cdh,
+                     double *sunRadVar_linke,
+                     double *sunRadVar_G_norm_extra,
+                     double *sunRadVar_alb);
 
 double brad_angle_loss(double, double *bh,
 		       struct SunGeometryVarDay *sunVarGeom,
@@ -49,13 +92,67 @@ void com_par(struct SunGeometryConstDay *sungeom,
 	     struct SunGeometryVarDay *sunVarGeom,
 	     struct GridGeometry *gridGeom,
 	     double latitude, double longitude);
+void com_par_parallel(double *sunGeom_lum_C11,
+                      double *sunGeom_lum_C13,
+                      double *sunGeom_lum_C22,
+                      double *sunGeom_lum_C31,
+                      double *sunGeom_lum_C33,
+                      double *sunGeom_sunrise_time,
+                      double *sunGeom_sunset_time,
+                      double *sunGeom_timeAngle,
+                      double *sunVarGeom_solarAltitude,
+                      double *sunVarGeom_sinSolarAltitude,
+                      double *sunVarGeom_tanSolarAltitude,
+                      double *sunVarGeom_solarAzimuth,
+                      double *sunVarGeom_sunAzimuthAngle,
+                      double *sunVarGeom_stepsinangle,
+                      double *sunVarGeom_stepcosangle,
+                      double *gridGeom_stepxy,
+                      double latitude);
+
 void com_par_const(double longitTime, struct SunGeometryConstDay *sungeom,
 		   struct GridGeometry *gridGeom);
+void com_par_const_parallel(double longitTime,
+                            double *sunGeom_lum_C11,
+                            double *sunGeom_lum_C13,
+                            double *sunGeom_lum_C22,
+                            double *sunGeom_lum_C31,
+                            double *sunGeom_lum_C33,
+                            double *sunGeom_sunrise_time,
+                            double *sunGeom_sunset_time,
+                            double *sunGeom_timeAngle,
+                            double *sunGeom_sindecl,
+                            double *sunGeom_cosdecl,
+                            double *gridGeom_sinlat,
+                            double *gridGeom_coslat);
+
 double lumcline2(struct SunGeometryConstDay *sungeom,
 		 struct SunGeometryVarDay *sunVarGeom,
 		 struct SunGeometryVarSlope *sunSlopeGeom,
 		 struct GridGeometry *gridGeom,
 		 unsigned char *horizonpointer);
+double lumcline2_parallel(double *sunGeom_timeAngle,
+                          int    *sunVarGeom_isShadow,
+                          double *sunVarGeom_solarAltitude,
+                          double *sunVarGeom_sunAzimuthAngle,
+                          double *sunVarGeom_z_orig,
+                          double *sunVarGeom_zmax,
+                          double *sunVarGeom_zp,
+                          double *sunVarGeom_tanSolarAltitude,
+                          double *sunVarGeom_stepsinangle,
+                          double *sunVarGeom_stepcosangle,
+                          double *sunSlopeGeom_longit_l,
+                          double *sunSlopeGeom_lum_C31_l,
+                          double *sunSlopeGeom_lum_C33_l,
+                          double *gridGeom_xx0,
+                          double *gridGeom_yy0,
+                          double *gridGeom_xg0,
+                          double *gridGeom_yg0,
+                          double *gridGeom_stepx,
+                          double *gridGeom_stepy,
+                          double *gridGeom_deltx,
+                          double *gridGeom_delty,
+                          unsigned char *horizonpointer);
 
 
 typedef double (*BeamRadFunc) (double sh, double *bh,

raster/r.sun/main.c

@@ -29,6 +29,11 @@ email: hofierka@geomodel.sk,marcel.suri@jrc.it,suri@geomodel.sk Thomas.Huld@jrc.
 /*v. 2.0 July 2002, NULL data handling, JH */
 /*v. 2.1 January 2003, code optimization by Thomas Huld, JH */
 /*v. 3.0 February 2006, several changes (shadowing algorithm, earth's curvature JH */
+/*v. 4.0 December 2016, OpenMP support Stanislav Zubal, Michal Lacko */
+
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
 
 #include <stdio.h>
 #include <stdlib.h>
@@ -43,6 +48,7 @@ email: hofierka@geomodel.sk,marcel.suri@jrc.it,suri@geomodel.sk Thomas.Huld@jrc.
 
 #include <grass/gis.h>
 #include <grass/raster.h>
+#include <grass/gmath.h>
 #include <grass/gprojects.h>
 #include <grass/glocale.h>
 #include "sunradstruct.h"
@@ -122,7 +128,11 @@ void joules2(struct SunGeometryConstDay *sunGeom,
 	     struct SolarRadVar *sunRadVar,
 	     struct GridGeometry *gridGeom,
 	     unsigned char *horizonpointer, double latitude,
-	     double longitude);
+	     double longitude,            
+	     double *Pbeam_e,
+	     double *Pdiff_e,
+	     double *Prefl_e,
+	     double *Pinsol_t);
 
 
 void calculate(double singleSlope, double singleAspect,
@@ -213,6 +223,8 @@ int main(int argc, char *argv[])
     double singleAlbedo;
     double singleLinke;
 
+    int threads;
+
     struct GModule *module;
     struct
     {
@@ -220,7 +232,7 @@ int main(int argc, char *argv[])
 	    *lin, *albedo, *longin, *alb, *latin, *coefbh, *coefdh,
 	    *incidout, *beam_rad, *insol_time, *diff_rad, *refl_rad,
 	    *glob_rad, *day, *step, *declin, *ltime, *dist, *horizon,
-	    *horizonstep, *numPartitions, *civilTime;
+	    *horizonstep, *numPartitions, *civilTime, *threads;
     }
     parm;
 
@@ -469,6 +481,15 @@ int main(int argc, char *argv[])
     parm.ltime->options = "0-24";
     parm.ltime->guisection = _("Time");
 
+    parm.threads = G_define_option();
+    parm.threads->key = "nprocs";
+    parm.threads->type = TYPE_INTEGER;
+    parm.threads->answer = "1";
+    parm.threads->options = "1-1000";
+    parm.threads->required = NO;
+    parm.threads->description =
+	_("Number of threads which will be used for parallel computing");
+
     /*
      * parm.startTime = G_define_option();
      * parm.startTime->key = "start_time";
@@ -551,6 +572,20 @@ int main(int argc, char *argv[])
 
     civiltime = parm.civilTime->answer;
 
+    sscanf(parm.threads->answer, "%d", &threads);
+    if (threads < 1)
+    {
+      G_warning(_("<%d> is not valid number of threads. Number of threads will be set to <%d>"),
+      threads, abs(threads));
+      threads = abs(threads);
+    }
+#if defined(_OPENMP)
+    omp_set_num_threads(threads);
+#else
+    threads = 1;
+#endif
+    G_message(_("Number of threads <%d>"), threads);
+
     if (civiltime != NULL) {
 	setUseCivilTime(TRUE);
 
@@ -790,7 +825,6 @@ int INPUT_part(int offset, double *zmax)
 {
     int finalRow, rowrevoffset;
     int numRows;
-    int numDigits;
     FCELL *cell1 = NULL, *cell2 = NULL;
     FCELL *cell3 = NULL, *cell4 = NULL, *cell5 = NULL, *cell6 = NULL, *cell7 =
 	NULL;
@@ -801,9 +835,8 @@ int INPUT_part(int offset, double *zmax)
 	fd7 = -1, row, row_rev;
     static int *fd_shad;
     int fr1 = -1, fr2 = -1;
-    int l, i, j;
+    int i, j;
     char *shad_filename;
-    char formatString[10];
     double angle_deg = 0.;
 
     finalRow = m - offset - m / numPartitions;
@@ -816,12 +849,10 @@ int INPUT_part(int offset, double *zmax)
     cell1 = Rast_allocate_f_buf();
 
     if (z == NULL) {
-	z = (float **)G_malloc(sizeof(float *) * (numRows));
-
-
-	for (l = 0; l < numRows; l++) {
-	    z[l] = (float *)G_malloc(sizeof(float) * (n));
-	}
+	if (!(z = G_alloc_fmatrix(numRows, n))) 
+        {
+            G_fatal_error(_("Out of memory"));
+        }
     }
 
     fd1 = Rast_open_old(elevin, "");
@@ -829,11 +860,10 @@ int INPUT_part(int offset, double *zmax)
     if (slopein != NULL) {
 	cell3 = Rast_allocate_f_buf();
 	if (s == NULL) {
-	    s = (float **)G_malloc(sizeof(float *) * (numRows));
-
-	    for (l = 0; l < numRows; l++) {
-		s[l] = (float *)G_malloc(sizeof(float) * (n));
-	    }
+            if (!(s = G_alloc_fmatrix(numRows, n))) 
+            {
+                G_fatal_error(_("Out of memory"));
+            }
 	}
 	fd3 = Rast_open_old(slopein, "");
 
@@ -843,11 +873,10 @@ int INPUT_part(int offset, double *zmax)
 	cell2 = Rast_allocate_f_buf();
 
 	if (o == NULL) {
-	    o = (float **)G_malloc(sizeof(float *) * (numRows));
-
-	    for (l = 0; l < numRows; l++) {
-		o[l] = (float *)G_malloc(sizeof(float) * (n));
-	    }
+	    if (!(o = G_alloc_fmatrix(numRows, n))) 
+            {
+                G_fatal_error(_("Out of memory"));
+            }
 	}
 	fd2 = Rast_open_old(aspin, "");
 
@@ -856,10 +885,10 @@ int INPUT_part(int offset, double *zmax)
     if (linkein != NULL) {
 	cell4 = Rast_allocate_f_buf();
 	if (li == NULL) {
-	    li = (float **)G_malloc(sizeof(float *) * (numRows));
-	    for (l = 0; l < numRows; l++)
-		li[l] = (float *)G_malloc(sizeof(float) * (n));
-
+            if (!(li = G_alloc_fmatrix(numRows, n))) 
+            {
+                G_fatal_error(_("Out of memory"));
+            }
 	}
 	fd4 = Rast_open_old(linkein, "");
     }
@@ -867,9 +896,10 @@ int INPUT_part(int offset, double *zmax)
     if (albedo != NULL) {
 	cell5 = Rast_allocate_f_buf();
 	if (a == NULL) {
-	    a = (float **)G_malloc(sizeof(float *) * (numRows));
-	    for (l = 0; l < numRows; l++)
-		a[l] = (float *)G_malloc(sizeof(float) * (n));
+            if (!(a = G_alloc_fmatrix(numRows, n))) 
+            {
+                G_fatal_error(_("Out of memory"));
+            }
 	}
 	fd5 = Rast_open_old(albedo, "");
     }
@@ -877,18 +907,22 @@ int INPUT_part(int offset, double *zmax)
     if (latin != NULL) {
 	cell6 = Rast_allocate_f_buf();
 	if (latitudeArray == NULL) {
-	    latitudeArray = (float **)G_malloc(sizeof(float *) * (numRows));
-	    for (l = 0; l < numRows; l++)
-		latitudeArray[l] = (float *)G_malloc(sizeof(float) * (n));
+            if (!(latitudeArray = G_alloc_fmatrix(numRows, n))) 
+            {
+                G_fatal_error(_("Out of memory"));
+            }
 	}
 	fd6 = Rast_open_old(latin, "");
     }
 
     if (longin != NULL) {
 	cell7 = Rast_allocate_f_buf();
-	longitudeArray = (float **)G_malloc(sizeof(float *) * (numRows));
-	for (l = 0; l < numRows; l++)
-	    longitudeArray[l] = (float *)G_malloc(sizeof(float) * (n));
+	if (latitudeArray == NULL) {
+            if (!(longitudeArray = G_alloc_fmatrix(numRows, n))) 
+            {
+                G_fatal_error(_("Out of memory"));
+            }
+	}
 
 	fd7 = Rast_open_old(longin, "");
     }
@@ -896,9 +930,10 @@ int INPUT_part(int offset, double *zmax)
     if (coefbh != NULL) {
 	rast1 = Rast_allocate_f_buf();
 	if (cbhr == NULL) {
-	    cbhr = (float **)G_malloc(sizeof(float *) * (numRows));
-	    for (l = 0; l < numRows; l++)
-		cbhr[l] = (float *)G_malloc(sizeof(float) * (n));
+            if (!(cbhr = G_alloc_fmatrix(numRows, n))) 
+            {
+                G_fatal_error(_("Out of memory"));
+            }
 	}
 	fr1 = Rast_open_old(coefbh, "");
     }
@@ -906,9 +941,10 @@ int INPUT_part(int offset, double *zmax)
     if (coefdh != NULL) {
 	rast2 = Rast_allocate_f_buf();
 	if (cdhr == NULL) {
-	    cdhr = (float **)G_malloc(sizeof(float *) * (numRows));
-	    for (l = 0; l < numRows; l++)
-		cdhr[l] = (float *)G_malloc(sizeof(float) * (n));
+            if (!(cdhr = G_alloc_fmatrix(numRows, n))) 
+            {
+                G_fatal_error(_("Out of memory"));
+            }
 	}
 	fr2 = Rast_open_old(coefdh, "");
     }
@@ -916,11 +952,11 @@ int INPUT_part(int offset, double *zmax)
     if (useHorizonData()) {
 	if (horizonarray == NULL) {
 	    horizonarray =
-		(unsigned char *)G_malloc(sizeof(char) * arrayNumInt *
-					  numRows * n);
+		(unsigned char *)G_calloc(arrayNumInt *
+					  numRows * n, sizeof(char));
 
-	    horizonbuf = (FCELL **) G_malloc(sizeof(FCELL *) * arrayNumInt);
-	    fd_shad = (int *)G_malloc(sizeof(int) * arrayNumInt);
+	    horizonbuf = (FCELL **) G_calloc(arrayNumInt, sizeof(FCELL *) );
+	    fd_shad = (int *)G_calloc(arrayNumInt, sizeof(int));
 	}
 	/*
 	 * if(ttime != NULL)
@@ -1298,7 +1334,11 @@ void joules2(struct SunGeometryConstDay *sunGeom,
 	     struct SunGeometryVarSlope *sunSlopeGeom,
 	     struct SolarRadVar *sunRadVar,
 	     struct GridGeometry *gridGeom, unsigned char *horizonpointer,
-	     double latitude, double longitude)
+	     double latitude, double longitude,
+	     double *Pbeam_e,
+	     double *Pdiff_e,
+	     double *Prefl_e,
+	     double *Pinsol_t)
 {
 
     double s0, dfr, dfr_rad;
@@ -1310,10 +1350,6 @@ void joules2(struct SunGeometryConstDay *sunGeom,
     double bh;
     double rr;
 
-    beam_e = 0.;
-    diff_e = 0.;
-    refl_e = 0.;
-    insol_t = 0.;
 
 
     com_par(sunGeom, sunVarGeom, gridGeom, latitude, longitude);
@@ -1326,21 +1362,21 @@ void joules2(struct SunGeometryConstDay *sunGeom,
 	if (sunVarGeom->solarAltitude > 0.) {
 	    if ((!sunVarGeom->isShadow) && (s0 > 0.)) {
 		ra = brad(s0, &bh, sunVarGeom, sunSlopeGeom, sunRadVar);	/* beam radiation */
-		beam_e += ra;
+		*Pbeam_e += ra;
 	    }
 	    else {
-		beam_e = 0.;
+		*Pbeam_e = 0.;
 		bh = 0.;
 	    }
 
 	    if ((diff_rad != NULL) || (glob_rad != NULL)) {
 		dra = drad(s0, bh, &rr, sunVarGeom, sunSlopeGeom, sunRadVar);	/* diffuse rad. */
-		diff_e += dra;
+		*Pdiff_e += dra;
 	    }
 	    if ((refl_rad != NULL) || (glob_rad != NULL)) {
 		if ((diff_rad == NULL) && (glob_rad == NULL))
 		    drad(s0, bh, &rr, sunVarGeom, sunSlopeGeom, sunRadVar);
-		refl_e += rr;	/* reflected rad. */
+		*Prefl_e += rr;	/* reflected rad. */
 	    }
 	}			/* solarAltitude */
     }
@@ -1377,9 +1413,9 @@ void joules2(struct SunGeometryConstDay *sunGeom,
 	    if (sunVarGeom->solarAltitude > 0.) {
 
 		if ((!sunVarGeom->isShadow) && (s0 > 0.)) {
-		    insol_t += dfr;
+		    *Pinsol_t += dfr;
 		    ra = brad(s0, &bh, sunVarGeom, sunSlopeGeom, sunRadVar);
-		    beam_e += dfr * ra;
+		    *Pbeam_e += dfr * ra;
 		    ra = 0.;
 		}
 		else {
@@ -1389,7 +1425,7 @@ void joules2(struct SunGeometryConstDay *sunGeom,
 		    dra =
 			drad(s0, bh, &rr, sunVarGeom, sunSlopeGeom,
 			     sunRadVar);
-		    diff_e += dfr * dra;
+		    *Pdiff_e += dfr * dra;
 		    dra = 0.;
 		}
 		if ((refl_rad != NULL) || (glob_rad != NULL)) {
@@ -1397,7 +1433,7 @@ void joules2(struct SunGeometryConstDay *sunGeom,
 			drad(s0, bh, &rr, sunVarGeom, sunSlopeGeom,
 			     sunRadVar);
 		    }
-		    refl_e += dfr * rr;
+		    *Prefl_e += dfr * rr;
 		    rr = 0.;
 		}
 	    }			/* illuminated */
@@ -1685,9 +1721,9 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
 
 
     if (incidout != NULL) {
-	lumcl = (float **)G_malloc(sizeof(float *) * (m));
+	lumcl = (float **)G_calloc((m), sizeof(float *));
 	for (l = 0; l < m; l++) {
-	    lumcl[l] = (float *)G_malloc(sizeof(float) * (n));
+	    lumcl[l] = (float *)G_calloc((n), sizeof(float *));
 	}
 	for (j = 0; j < m; j++) {
 	    for (i = 0; i < n; i++)
@@ -1696,9 +1732,9 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
     }
 
     if (beam_rad != NULL) {
-	beam = (float **)G_malloc(sizeof(float *) * (m));
+	beam = (float **)G_calloc((m), sizeof(float *));
 	for (l = 0; l < m; l++) {
-	    beam[l] = (float *)G_malloc(sizeof(float) * (n));
+	    beam[l] = (float *)G_calloc((n), sizeof(float *));
 	}
 
 	for (j = 0; j < m; j++) {
@@ -1708,9 +1744,9 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
     }
 
     if (insol_time != NULL) {
-	insol = (float **)G_malloc(sizeof(float *) * (m));
+	insol = (float **)G_calloc((m), sizeof(float *));
 	for (l = 0; l < m; l++) {
-	    insol[l] = (float *)G_malloc(sizeof(float) * (n));
+	    insol[l] = (float *)G_calloc((n), sizeof(float *));
 	}
 
 	for (j = 0; j < m; j++) {
@@ -1720,9 +1756,9 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
     }
 
     if (diff_rad != NULL) {
-	diff = (float **)G_malloc(sizeof(float *) * (m));
+	diff = (float **)G_calloc((m), sizeof(float *));
 	for (l = 0; l < m; l++) {
-	    diff[l] = (float *)G_malloc(sizeof(float) * (n));
+	    diff[l] = (float *)G_calloc((n), sizeof(float *));
 	}
 
 	for (j = 0; j < m; j++) {
@@ -1732,9 +1768,9 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
     }
 
     if (refl_rad != NULL) {
-	refl = (float **)G_malloc(sizeof(float *) * (m));
+	refl = (float **)G_calloc((m), sizeof(float *));
 	for (l = 0; l < m; l++) {
-	    refl[l] = (float *)G_malloc(sizeof(float) * (n));
+	    refl[l] = (float *)G_calloc((n), sizeof(float *));
 	}
 
 	for (j = 0; j < m; j++) {
@@ -1744,9 +1780,9 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
     }
 
     if (glob_rad != NULL) {
-	globrad = (float **)G_malloc(sizeof(float *) * (m));
+	globrad = (float **)G_calloc((m), sizeof(float *));
 	for (l = 0; l < m; l++) {
-	    globrad[l] = (float *)G_malloc(sizeof(float) * (n));
+	    globrad[l] = (float *)G_calloc((n), sizeof(float *));
 	}
 
 	for (j = 0; j < m; j++) {
@@ -1777,7 +1813,7 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
     else {
 	setTimeOffset(0.);
     }
-
+    int shadowoffset_base = shadowoffset;
     for (j = 0; j < m; j++) {
 	G_percent(j, m - 1, 2);
 
@@ -1788,10 +1824,12 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
 
 	}
 	sunVarGeom.zmax = zmax;
-
-
+        shadowoffset_base = (j % (numRows)) * n * arrayNumInt;
+    #pragma omp parallel firstprivate(q1,tan_lam_l,z1,i,shadowoffset,longitTime,coslat,coslatsq,func,latitude,longitude,sin_phi_l,latid_l,sin_u,cos_u,sin_v,cos_v,lum, gridGeom,sunGeom,sunVarGeom,sunSlopeGeom,sunRadVar, elevin,aspin,slopein,civiltime,linkein,albedo,latin,coefbh,coefdh,incidout,longin,horizon,beam_rad,insol_time,diff_rad,refl_rad,glob_rad,mapset,per,decimals,str_step )
+    {
+      #pragma omp for schedule(dynamic) 
 	for (i = 0; i < n; i++) {
-
+            shadowoffset = shadowoffset_base + (arrayNumInt * i);
 	    if (useCivilTime()) {
 		/* sun travels 15deg per hour, so 1 TZ every 15 deg and 15 TZs * 24hrs = 360deg */
 		longitTime = -longitudeArray[arrayOffset][i] / 15.;
@@ -1922,27 +1960,31 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
 		    }
 		    else lumcl[j][i] = UNDEFZ;
 		}
+ 
+ 		double Pbeam_e = 0.;
+ 		double Pdiff_e = 0.;
+ 		double Prefl_e = 0.;
+ 		double Pinsol_t = 0.;
 
 		if (someRadiation) {
 		    joules2(&sunGeom, &sunVarGeom, &sunSlopeGeom, &sunRadVar,
 			    &gridGeom, horizonarray + shadowoffset, latitude,
-			    longitude);
+  			    longitude, &Pbeam_e, &Pdiff_e, &Prefl_e, &Pinsol_t);
 		    if (beam_rad != NULL)
-			beam[j][i] = (float)beam_e;
+ 			beam[j][i] = (float)Pbeam_e;
 		    if (insol_time != NULL)
-			insol[j][i] = (float)insol_t;
+ 			insol[j][i] = (float)Pinsol_t;
 		    /*  G_debug(3,"\n %f",insol[j][i]); */
 		    if (diff_rad != NULL)
-			diff[j][i] = (float)diff_e;
+ 			diff[j][i] = (float)Pdiff_e;
 		    if (refl_rad != NULL)
-			refl[j][i] = (float)refl_e;
+ 			refl[j][i] = (float)Prefl_e;
 		    if (glob_rad != NULL)
-			globrad[j][i] = (float)(beam_e + diff_e + refl_e);
+ 			globrad[j][i] = (float)(Pbeam_e + Pdiff_e + Prefl_e);
 		}
 
 	    }			/* undefs */
-	    shadowoffset += arrayNumInt;
-	}
+	}}
 	arrayOffset++;
     }
 

raster/r.sun/testsuite/test_rsun.py

@@ -0,0 +1,117 @@
+# -*- coding: utf-8 -*-
+from grass.gunittest.case import TestCase
+from grass.gunittest.main import test
+from grass.gunittest.gmodules import SimpleModule
+
+
+class TestRSunMode2(TestCase):
+    elevation = 'elevation'
+    elevation_attrib = 'elevation_attrib'
+    elevation_threads = 'elevation_threads'
+    slope = 'rsun_slope'
+    aspect = 'rsun_aspect'
+    beam_rad = 'beam_rad'
+    glob_rad = 'glob_rad'
+    insol_time = 'insol_time'
+    glob_rad_threads = 'glob_rad_threads'
+
+    @classmethod
+    def setUpClass(cls):
+        cls.use_temp_region()
+        cls.runModule('g.region', n=223500, s=220000, e=640000, w=635000, res=10)
+        cls.runModule('r.slope.aspect', elevation=cls.elevation, slope=cls.slope, aspect=cls.aspect)
+
+    @classmethod
+    def tearDownClass(cls):
+        cls.del_temp_region()
+        cls.runModule('g.remove', type=['raster'],
+                      name=[cls.slope, cls.aspect, cls.insol_time, cls.beam_rad, cls.glob_rad, cls.glob_rad_threads], flags='f')
+
+    def setUp(self):
+        self.rsun = SimpleModule('r.sun', elevation=self.elevation, slope=self.slope, aspect=self.aspect,
+                                 day=172, beam_rad=self.beam_rad, glob_rad=self.glob_rad, insol_time=self.insol_time,
+                                 overwrite=True)
+
+    def test_more_threads(self):
+        self.assertModule(self.rsun)
+        try:
+            self.rsun.inputs['nprocs'].value = 4
+            self.rsun.outputs.glob_rad = self.glob_rad_threads
+            self.assertModule(self.rsun)
+            self.assertRastersNoDifference(self.glob_rad, self.glob_rad_threads, precision=1e-8)
+        except KeyError:
+            # original version of r.sun without parallel processing
+            return
+
+    def test_run_outputs(self):
+        self.assertModule(self.rsun)
+        self.assertRasterExists(name=self.beam_rad)
+        self.assertRasterExists(name=self.glob_rad)
+        self.assertRasterExists(name=self.insol_time)
+        # beam_rad
+        values = 'min=6561.3505859375\nmax=7680.93505859375\nmean=7638.91454059886'
+        self.assertRasterFitsUnivar(raster=self.beam_rad, reference=values, precision=1e-8)
+        # glob_rad
+        values = 'min=7829.4921875\nmax=8889.4765625\nmean=8851.4872842213'
+        self.assertRasterFitsUnivar(raster=self.glob_rad, reference=values, precision=1e-8)
+        # insol_time
+        values = 'min=11.5\nmax=14\nmean=13.8521038175691'
+        self.assertRasterFitsUnivar(raster=self.insol_time, reference=values, precision=1e-8)
+
+
+class TestRSunMode1(TestCase):
+    elevation = 'elevation'
+    elevation_attrib = 'elevation_attrib'
+    elevation_threads = 'elevation_threads'
+    slope = 'rsun_slope'
+    aspect = 'rsun_aspect'
+    beam_rad = 'beam_rad'
+    glob_rad = 'glob_rad'
+    incidout = 'incidout'
+    glob_rad_threads = 'glob_rad_threads'
+
+    @classmethod
+    def setUpClass(cls):
+        cls.use_temp_region()
+        cls.runModule('g.region', n=223500, s=220000, e=640000, w=635000, res=10)
+        cls.runModule('r.slope.aspect', elevation=cls.elevation, slope=cls.slope, aspect=cls.aspect)
+
+    @classmethod
+    def tearDownClass(cls):
+        cls.del_temp_region()
+        cls.runModule('g.remove', type=['raster'],
+                      name=[cls.slope, cls.aspect, cls.incidout, cls.beam_rad, cls.glob_rad, cls.glob_rad_threads], flags='f')
+
+    def setUp(self):
+        self.rsun = SimpleModule('r.sun', elevation=self.elevation, slope=self.slope, aspect=self.aspect,
+                                 day=172, time=18, beam_rad=self.beam_rad, glob_rad=self.glob_rad, incidout=self.incidout,
+                                 overwrite=True)
+
+    def test_more_threads(self):
+        self.assertModule(self.rsun)
+        try:
+            self.rsun.inputs['nprocs'].value = 4
+            self.rsun.outputs.glob_rad = self.glob_rad_threads
+            self.assertModule(self.rsun)
+            self.assertRastersNoDifference(self.glob_rad, self.glob_rad_threads, precision=1e-8)
+        except KeyError:
+            # original version of r.sun without parallel processing
+            return
+
+    def test_run_outputs(self):
+        self.assertModule(self.rsun)
+        self.assertRasterExists(name=self.beam_rad)
+        self.assertRasterExists(name=self.glob_rad)
+        self.assertRasterExists(name=self.incidout)
+        # beam_rad
+        values = 'min=0\nmax=355.780944824219\nmean=124.615522080923'
+        self.assertRasterFitsUnivar(raster=self.beam_rad, reference=values, precision=1e-8)
+        # glob_rad
+        values = 'min=32.1095008850098\nmax=451.527923583984\nmean=178.057346498427'
+        self.assertRasterFitsUnivar(raster=self.glob_rad, reference=values, precision=1e-8)
+        # incidout
+        values = 'min=0.0100772567093372\nmax=40.5435371398926\nmean=13.2855086254291'
+        self.assertRasterFitsUnivar(raster=self.incidout, reference=values, precision=1e-8)
+
+if __name__ == '__main__':
+    test()