Removing double*

DESCRIPTION

@@ -2,7 +2,7 @@ Package: medfate
 Type: Package
 Title: Mediterranean Forest Simulation
 Version: 4.8.0
-Date: 2024-10-24
+Date: 2024-10-17
 Authors@R: c(
 	person('Miquel', 'De Cáceres', role=c('aut','cre', 'cph'),
 	email='miquelcaceres@gmail.com', comment = c(ORCID = "0000-0001-7132-2080")),

NEWS.md

@@ -1,7 +1,6 @@
 # medfate 4.8.0
 * Code optimization: LAI distribution
-* Using ANSI C in internal vectors
-* Global communication structures
+* Communication structures to improve memory usage in medfateland
 
 # medfate 4.7.0
 * LAI can be estimated without the competition effect of cohorts of larger size

src/hydraulics.cpp

@@ -234,17 +234,15 @@ double Egamma(double psi, double kxylemmax, double c, double d, double psiCav =
   else if(psi==0.0) return(0.0);
   double h = 1.0/c;
   double z = pow(psi/d,c);
-  double* pq = incgam(h,z);
+  NumericVector pq = incgam(h,z);
   double g = tgamma(h)*pq[0]; //Upper incomplete gamma, without the normalizing factor
-  delete[] pq;
   double E = kxylemmax*(-d/c)*g;
   if(psiCav<0.0) { //Decrease E from 0 to psiCav (avoid recursiveness!)
     if(psiCav < psi) {
       E = xylemConductance(psiCav,kxylemmax,c,d)*(-psi); //square integral
     } else {
-      double* pq = incgam(h,pow(psiCav/d,c));
+      NumericVector pq = incgam(h,pow(psiCav/d,c));
       double Epsimin = kxylemmax*(-d/c)*tgamma(h)*pq[0];
-      delete[] pq;
       E = E - Epsimin + xylemConductance(psiCav,kxylemmax,c,d)*(-psiCav); //Remove part of the integral corresponding to psimin and add square integral
     }
   }

src/incgamma.cpp

@@ -171,7 +171,7 @@ double lnec(double x) {
 //   chepolsum=a(0)/2.0_r8-r+h*x
 //   ENDIF
 //   END FUNCTION chepolsum
-double chepolsum(double x, double* a, int n) {
+double chepolsum(double x, NumericVector a, int n) {
   if(n==0) {
     return(a[0]/2.0);
   } else if (n==1) {
@@ -220,7 +220,7 @@ double chepolsum(double x, double* a, int n) {
 // ENDIF
 //   END FUNCTION auxgam
 double auxgam(double x) {
-  double dr[18];
+  NumericVector dr(18);
   double auxgamm;
   if(x<0.0) {
     auxgamm = -(1.0+(1.0+x)*(1.0+x)*auxgam(1.0+x))/(1.0-x);
@@ -316,8 +316,8 @@ double lngam1(double x) {
 //   END FUNCTION stirling
 double  stirling(double x) {
   double stirling, z;
-  double a[18];
-  double c[7];
+  NumericVector a(18);
+  NumericVector c(7);
   if(x<dwarf) {
     stirling = giant; 
   } else if(x<1.0) {
@@ -1092,7 +1092,7 @@ double qfraction(double a, double x, double dp){
 //     ENDIF
 //   ENDIF
 //   END SUBROUTINE incgam
-double* incgam(double a, double x) {
+NumericVector incgam(double a, double x) {
   double lnx, p = NA_REAL, q = NA_REAL;
   double dp;
   if(x<dwarf) {
@@ -1139,10 +1139,7 @@ double* incgam(double a, double x) {
       }
     }
   }
-  double* res = new double[2];
-  res[0] = p;
-  res[1] = q;
-  return(res);
+  return(NumericVector::create(p,q));
 }
 
 
@@ -1831,17 +1828,15 @@ double invincgam(double a, double p, double q) {
       } else {
         r=exp(dlnr);
         if(pcase) {
-          double* pq = incgam(a,x);
+          NumericVector pq = incgam(a,x);
           px = pq[0];
           qx = pq[1];
           ck[0]=-r*(px-p);
-          delete[] pq;
         } else {
-          double* pq = incgam(a,x);
+          NumericVector pq = incgam(a,x);
           px = pq[0];
           qx = pq[1];
           ck[0]=r*(qx-q);
-          delete[] pq;
         }
         ck[1]=(x-a+1.0)/(2.0*x);
         ck[2]=(2.0*x2-4.0*x*a+4.0*x+2.0*a2-3.0*a+1.0)/(6.0*x2);
@@ -1862,17 +1857,15 @@ double invincgam(double a, double p, double q) {
       fp=-sqrt(a/twopi)*exp(-0.5*a*y*y)/(gamstar(a));
       r=-(1.0/fp)*x;
       if(pcase) {
-        double* pq = incgam(a,x);
+        NumericVector pq = incgam(a,x);
         px = pq[0];
         qx = pq[1];
         ck[0]=-r*(px-p);
-        delete[] pq;
       } else {
-        double* pq = incgam(a,x);
+        NumericVector pq = incgam(a,x);
         px = pq[0];
         qx = pq[1];
         ck[0]=r*(qx-q);
-        delete[] pq;
       }
       ck[1]=(x-a+1.0)/(2.0*x);
       ck[2]=(2.0*x2-4.0*x*a+4.0*x+2.0*a2-3.0*a+1.0)/(6.0*x2);

src/incgamma.h

@@ -5,6 +5,6 @@
 #endif
 using namespace Rcpp;
 
-double* incgam(double a, double x);
+NumericVector incgam(double a, double x);
 double invincgam(double a, double p, double q);
 double errorfunction(double x, bool erfcc, bool expo);

src/soil_thermodynamics.cpp

@@ -188,14 +188,14 @@ NumericVector temperatureChange(NumericVector widths, NumericVector Temp,
   NumericVector a(nlayers, 0.0), b(nlayers, 0.0), c(nlayers, 0.0), d(nlayers, 0.0), e(nlayers, 0.0), f(nlayers, 0.0);
 
   //Estimate layer interfaces
-  double* dZ_m = new double[nlayers];
+  NumericVector dZ_m(nlayers);
   for(int l=0;l<nlayers;l++) dZ_m[l] = widths[l]*0.001; //mm to m
 
-  double* dZUp = new double[nlayers];
-  double* dZDown = new double[nlayers];
-  double* Zcent = new double[nlayers];
-  double* Zup = new double[nlayers];
-  double* Zdown = new double[nlayers];
+  NumericVector dZUp(nlayers,0.0);
+  NumericVector dZDown(nlayers,0.0);
+  NumericVector Zcent(nlayers,0.0);
+  NumericVector Zup(nlayers,0.0);
+  NumericVector Zdown(nlayers,0.0);
   for(int l=0;l<nlayers;l++) {
     if(l==0) { //first layer
       dZUp[l] = dZ_m[0]/2.0; //Distance from ground to mid-layer
@@ -214,8 +214,8 @@ NumericVector temperatureChange(NumericVector widths, NumericVector Temp,
       dZDown[l] = dZ_m[l]/2.0;
     }
   }
-  double* k_up = new double[nlayers];
-  double* k_down = new double[nlayers];
+  NumericVector k_up(nlayers,0.0);
+  NumericVector k_down(nlayers,0.0);
   for(int l=0;l<nlayers;l++) {
     k_up[l] = 0.0;
     k_down[l] = 0.0;