/* Program weevils.c */ /* LAST MODIFIED 6 Sept, 2007 */ /* This version of the code uses integer numbers of weevils, egg holes, * etc., rather than real values based on averages (expected values). * This should be more accurate, since these numbers typically are * small. Also, this approach allows weevil numbers to be zero, not * possible re expected values. As in the Nov, 05 code from which * this code is derived, this code generates starting seeds randomly * rather than reading them from the parameter list for the sims. * * This is program weevils.c, a 'C' program simulating * weevil-spruce interaction when using TWO, SIX, EIGHTEEN or THIRTY * clones. With minor modifications, this code can be used when randomly * generating ANY number of clones. * First draft by Kyle Schertzer, Seth McElvaney and John Bishir, 1997-1998 * 1998 may 5 - Reformatted by Ken Polsson, removing C++ extensions, * preparing for use in TASS. * * 16 December, 1998 Modified by Ken Polsson and John Bishir * and incorporated as a subroutine in TASS * * April, 2002, onward -- Further modifications by Ken Polsson and John Bishir * * November, 2004 Modifications by Ken Polsson and John Bishir * * April, 2005 Further modifications by Ken Polsson and John Bishir * * August 14, 2006 Added chance_check() * September 12, 2006 Added PoissonRandom() */ #include #include #include #include #include "random.h" #include "weevils.h" #include "define.h" #include "plot.h" #include "species.h" #include "tree.h" #include "merch.h" /* #define DEBUG 1 */ #define FALSE 0 #define TRUE 1 #define izero 0 /* QUANTITIES ASSOCIATED WITH THE PLANTATION AS A WHOLE */ static int noweevil; /* set noweevil equal 1 to have no weevils in the stand and equal 0 to have weevils*/ static int numcln; /* number of clones */ static int nr; /* total number of rows, with nr = nrblk*rpb */ static int nc; /* total number of columns, with nc = ncblk*cpb */ /* Note -- 'Edge' effects are minimized by choosing blocks as nearly square as possible */ #define GRID 3.0 /* grid distance between trees (in meters - used to * transform wind speed in km/hr into grid lines per second) */ /* * Values of WIND PARAMETERS * * In the simulation the wind vector will be randomly generated each day. * Speed is a gamma(WLAMBDA, WTEE)random variable; wind speed is first * specified in kilometers per hour, and then transformed into number of * grid lines per second. For now, though, I will use a uniform distribution * of wind speed */ #define WMEAN 5.0 /* mean wind speed (km/hr) (when I use a gamma dist) */ #define WVAR 25.0 /* variance of wind speed (re gamma dist) */ #define WLAMBDA 0.2 /* first parameter in a gamma distribution, * set it to WMEAN/WVAR */ #define WTEE 1.0 /* second param in a gamma distribution; set it to * WLAMBDA*WMEAN */ #define WNS 0.75 /* proportion of time wind direction is from south (in a * negative direction) */ #define WEW 0.75 /* proportion of time wind direction is from west */ static float urand; /* Uniform random deviate */ /* CONSTANTS AND VARIABLES ASSOCIATED WITH TEMPERATURE */ /* FIXED QUANTITIES */ #define TVAR 12.25 /* variance of temp around its mean, tavg */ #define TDVAR 16.0 /* var of tdev around its mean, tdmean */ #define THRESH 7.2 /* temperature threshold below which no weevil activity * occurs, and no degree-days are accumulated */ #define DEG1 0.002 /* constant governing winter survival rate */ static float tmean; /* yearly mean temperature (oC) */ #define TAMPL 7.25 /* amplitude of maximum deviation of average temp from * yearly mean; computed as July average minus tmean */ /* QUANTITIES PRODUCED BY TASS; NEW VALUES EACH YEAR */ static float **height; /* tree height */ static float **crown; /* crown volume */ static float **length; /* length of leader (cm) */ static float **diam; /* diameter of last year's leader */ /* CONSTANTS THAT ARE RANDOMLY GENERATED AT TIME OF 'PLANTING' */ static float **r0; /* resin flow rate (local - not related to * total resin in the tree) */ static float **a0; /* attractiveness of odors produced in the * crown and through volatilization of resin */ static float **t0; /* toxicity of resin */ static float **g0; /* * PARAMETERS USED IN GENERATING r0, t0, and a0 * All three are normally distributed; we only need relative values, so * units are chosen so these quantities are on comparable scales */ #define A0MEAN 5.0 /* mean of distribution for a0 */ #define A0VAR 0.25 /* variance of distribution for a0 */ #define R0MEAN 5.0 /* mean of normal distribution for generating r0 */ #define R0VAR 2.25 /* variance of distribution for r0 */ #define T0MEAN 5.0 /* mean of distribution for t0 */ #define T0VAR 2.25 /* variance of distribution for t0 */ /* PARAMETERS USED IN THE MODEL EQUATIONS */ /* Except for volatiles, rates are per day, in the summer */ /* PARAMETERS IN THE EQUATION FOR ADULTS (NONOVIPOSITORS) */ #define DADULT 0.02 /* adult mortality rate */ #define RAD1 0.015 /* governs rate at which non-ovipositors are * repelled their present tree, in absence of * both resin flow and resin toxicity */ #define RAD2 0.1 /* governs change in rate at which non-ovipositors are * repelled from their present tree, per unit increase * in resin flow rate */ #define RAD3 0.1 /* governs change in rate at which non-ovipositors are * repelled from their present tree, per unit increase * in resin toxicity */ #define RAD4 60.0 /* governs change in rate at which non-ovipositors are * repelled from their present tree, as a function of * number of holes */ #define RADDIFFUS 0.16 /* adult (non-ovipositor) diffusion rate */ static float wstay; /* the probability a weevil stays on a tree and does * not diffuse randomly; the probability of diffusing * randomly is (1-wstay) */ /* PARAMETERS IN THE EQUATION FOR OVIPOSITORS */ #define DOVI 0.02 /* mortality rate for ovipositing weevils */ #define ROV1 0.07 /* governs rate at which ovipositors cease laying eggs, in * absence of resin flow, toxicity, and existing holes */ #define ROV2 0.25 /* governs relative rate at which ovipositors cease laying * eggs, per unit increase in resin flow rate */ #define ROV3 0.25 /* same as ROV2, but related to resin toxicity */ #define ROV4 1.0 /* same as ROV2, but re number of holes */ #define RNOV1 0.5 /* governs rate at which non-ovipositors begin laying * eggs, in absence of resin flow, toxicity, * and existing holes */ #define RNOV2 0.01 /* governs relative decrease in rate at which * non-ovipositors begin laying eggs, per * unit increase in resin flow rate */ #define RNOV3 0.01 /* same as RNOV2, but related to resin toxicity */ #define RNOV4 1.0 /* same as RNOV2, but re number of holes */ #define multre1 14.7 /* multiplier for the rate of oviposition per ovipositor; * used to calibrate amount of leader damage caused by * weevils and, indirectly, the amount of merchantable * timber produced in a stand */ /* PARAMETERS IN THE EQUATION FOR HOLES */ #define RH1 2.0 /* rate at which non-ovipositors create feeding holes, in * absence of resin, toxicity, or other holes */ #define RH3 0.06 /* proportionality constant - relates to the relative * decrease in rate of hole creation, per unit * increase in r0 */ #define RH4 0.06 /* same as RH3, but re toxicity instead of flow */ #define RH5 0.14 /* relative rate of attack on a leader when there are holes * already in a leader of known length and diameter */ /* PARAMETERS IN THE EQUATION FOR CONCENTRATION OF VOLATILES */ #define RVDECOMP 0.06 /* volatile decomposition rate (includes vertical * diffusion) */ #define RVDIFFUS 3.2 /* volatile horizontal diffusion rate */ #define RVCREATE 1.0 /* volatile creation rate, per unit crown volume * The value of parameter RVCREATE is arbitrary; * it appears only in ratios and thus does * not affect the computed values */ static float vstay; /* the proportion of volatiles that stay 'on' a tree and do * not diffuse randomly; the proportion diffusing randomly * is (1-vstay) */ /* Variables used in i,j loops for odors and volatiles */ static float **int1v; static float **int2v; static float **int3v; static float tol; /* VARIABLES THAT INDICATE THE STATUS OF TREES AND THEIR WEEVIL POPULATIONS */ /* updated each day -- "new" = today, "old" = yesterday * The i,j entry in any matrix indicates the magnitude - on tree i,j - of the variable associated with that matrix */ static int **livebrood; /* Number of surviving brood-of-the-year adults */ static int **newad; /* Number of adult weevils, not including brood adults, but including ovipositors */ static int **oldad; static int **newovi; /* Number of ovipositing weevils */ static int **oldovi; static int ***newjuv; /* Number of juvenile weevils */ static int ***oldjuv; static int **maxjuv; /* Largest number of juvenile weevils */ static int **emerge; /* Total number of mature juveniles that emerge * on a particular day */ static float ***newdeg; /* Number of degree-days accumulated */ static float ***olddeg; static int **newh; /* Total number of weevil holes in leader */ static int **oldh; static int **neweh; /* Number of egg holes */ static int **oldeh; static float **newvol; /* Volatile concentration */ static float **oldvol; static float **create; /* Amount of volatiles created in crown */ static int **repel; /* Number of adult weevils repelled today */ static int **postlarv2; /* Number of juveniles past larval stage 2 in each tree */ static int **youngad; /* Number of juveniles past pupal stage in each tree */ static int **numegg; /* Number of eggs in each tree */ static int **totnewegg; /* Total # of eggs laid in each tree this year */ static int **totemerge; /* Total # of brood adults from each tree, this year */ static int **numlarv; /* Number of larvae in each tree */ static int **numpupa; /* Number of pupae in each tree */ static int **totjuv; /* Current total number of juveniles in each tree */ static int **randif; /* Number of randomly diffusing weevils */ static float **damage; /* Weevil damage to leaders this year */ static int **leader; /* leader = 1 if leader is alive; otherwise, leader = 0 */ /* MISCELLANEOUS OTHER VARIABLES */ static float xxx; static int rancln; /* Variable 'rancln' indicates whether clone parameters * a0, r0, t0, and g0 are specified in the code * (indicated by rancln = 0) or are generated * randomly (indicated by rancln = 1 */ static int nsim; static int pattern; /* Variable 'pattern' indicates which type of planting pattern is used. pattern = 0 indicates a pattern in which each tree is a separately generated clone. There are as many clones as trees (nr x nc clones altogether) pattern = 11 indicates a pattern in which seedlings are planted in a random mixture of ramets of a number (numcln) of already generated clones. This pattern is always used when there is only one clone (numcln = 1). pattern = 12 indicates a pattern in which trees are planted in a mosaic of clonal blocks, using a random arrangement of already generated clones. Each clone appears in the same number (bpc) of blocks. pattern = 21 indicates a planting pattern consisting of two clonal blocks, one resistant and one susceptible pattern = 22 indicates a planting pattern involving two clones, one resistant and one susceptible, in a 2x2 pattern of blocks pattern = 23 indicates a planting pattern involving two clones, one resistant and one susceptible, in a 4x4 pattern of blocks pattern = 61 indicates a planting pattern involving 6 clones, one in each of 6 blocks pattern = 62 indicates a planting pattern involving 6 clones, in a 6x6 pattern of blocks; details are in the code below pattern = 181 indicates a planting pattern involving 18 clones, one in of 18 blocks pattern = 183 indicates a planting pattern involving 18 clones, each clone being used in 3 randomly placed blocks pattern = 301 indicates a planting pattern involving 30 clones, one in of 30 blocks pattern = 303 indicates a planting pattern involving 30 clones, each clone being used in 3 randomly placed blocks */ static int *clntree; /* clntree[m] = yearly # of live trees of clone type m */ static int *startclntree; /* startclntree[m] = # of stems planted of clone type m */ static int **wee_tree; /* wee_tree[i][j] = 1 if the tree in position ij * is alive, and = 0, if dead */ static float **wmerch; /* wmerch[i][j] = merchantable volume in tree ij */ static float summerch; /* total merchantable volume in the stand */ static int **adjacent; /* adjacent[i][j] = number of live trees * adjacent to tree i,j */ static int **treetype; /* Re CASE1A, treetype[i][j] indicates * the number of the clone in position i,j */ static float *startw; /* startw[i] = average number of weevils per tree in trees *of clone type i at the start of the season */ static int weedist; /* Indicates the initial distribution of weevils in the stand weedist = 1 indicates INITDIST1 in the code, in which the initial number of weevils is the same on each tree weedist = 2 indicates INITDIST2 in the code, in which 'wee2' weevils are placed on each of 'corner' trees in the upper left corner of a stand, with no weevils on other trees, the numbers numwee and corner being chosen by the user */ static float approxha; /* approximate area of the stand */ static int repnum; /* the number of the rep, in repeated runs */ static float g0mean; /* mean of g0 - this is the average site index for * this species in this stand location */ /* The next 4 items serve as multipliers to check on respective sensitivities of volumme to the parameters a0, g0, r0 and t0 */ static float a0mult; static float g0mult; static float r0mult; static float t0mult; /* The next four quantities provide temporary storage for clone parameters * (planting); 50 is maximum number of clones to be used */ static float t0cl[51]; /* t0cl[m] = resin toxicity for clone # m */ static float a0cl[51]; /* a0cl[m] = crown odor attractiveness of clone # m */ static float r0cl[51]; /* r0cl[m] = resin flow rate for clone # m */ static float g0cl[51]; /* TASS tree height vigor - the average site index for clone m * The actual site index for a ramet depends on its site */ static long int firstseed = 0; /* used for setting up random number streams */ static long int secondseed = 0; struct weevil_parameters *weevil_move; static short chance_check (short, float); static void weevil_plant (FILE *, struct weevil_parameters *); static void weevil_plant (FILE *errfile, struct weevil_parameters *params) { /* Parameters re tree growth rate; g0 (tree vigour) and r0 (resin flow rate) are correlated */ float g0sd; /* standard deviation of g0 */ float g0var; /* variance of g0 */ float rho; /* coeff of correlation between r0 and g0 */ float r0sd; /* std dev of r0 */ float g0r0var; /* conditional variance of g0, given r0 */ float g0r0mean; /* conditional mean of g0, given r0 */ float sitevar; /* Site variance */ float hsq; /* Heritibility - the ratio g0var/(g0var + sitevar) */ int clonenumber[100]; /* Serial list of clones to be placed in the blocks */ int **blocktype; /* Re CASE1B, blocktype[i][j] is the number of the clone * in the i-th row of blocks and j-th column of blocks */ int *clonelist; /* Random arrangement of chosen clones */ int randcl; /* randcl = 1 if clones are chosen randomly and randcl = 0 if * clone params are specified by the user */ int bpc; /* number of blocks per clone. When using Planting Pattern 11 * or Planting Pattern 0, bpc should be set to 0 */ int nrblk; /* number of rows of clone blocks */ int ncblk; /* number of columns of clone blocks */ /* Must have numcln*bpc = nrblk*ncblk */ int rpb; /* number of rows per clone block */ int cpb; /* number of columns per clone block */ int num1, num2, num3, num4 ; int itn; /* GENERAL PURPOSE LOOP VARIABLES */ int i, j, k, m, n ; #ifdef DEBUG printf ("weevils.c: weevil_plant(): begin\n"); #endif /* get parameters */ nsim = params->nsim; noweevil = params->noweevil; pattern = params->pattern; tmean = params->tmean; weedist = params->weedist; numcln = params->numcln; nrblk = params->nrblk; ncblk = params->ncblk; rpb = params->rpb; cpb = params->cpb; bpc = params->bpc; rancln = params->rancln; repnum = params->repnum; approxha = params->approxha; a0mult = params->a0mult; g0mult = params->g0mult; r0mult = params->r0mult; t0mult = params->t0mult; /* Replace these two lines for obtaining initial seeds with the new * lines below */ firstseed = params->rand_seed1; secondseed = params->rand_seed2; /* Generate initial seed values randomly */ /* srand(time(NULL)); firstseed = rand(); srand(time(NULL)); secondseed = rand(); */ nr = nrblk * rpb; nc = ncblk * cpb; /* LOTS OF ARRAYS TO INITIALIZE */ crown = (float **) malloc ((nr+2) * sizeof (float *)); for (i = 0; i < nr+2; i++) crown[i] = (float *) malloc ((nc+2) * sizeof (float)); length = (float **) malloc ((nr+2) * sizeof (float *)); for (i = 0; i < nr+2; i++) length[i] = (float *) malloc ((nc+2) * sizeof (float)); height = (float **) malloc ((nr+2) * sizeof (float *)); for (i = 0; i < nr+2; i++) height[i] = (float *) malloc ((nc+2) * sizeof (float)); diam = (float **) malloc ((nr+2) * sizeof (float *)); for (i = 0; i < nr+2; i++) diam[i] = (float *) malloc ((nc+2) * sizeof (float)); r0 = (float **) malloc ((nr+2) * sizeof (float *)); for (i = 0; i < nr+2; i++) r0[i] = (float *) malloc ((nc+2) * sizeof (float)); a0 = (float **) malloc ((nr+2) * sizeof (float *)); for (i = 0; i < nr+2; i++) a0[i] = (float *) malloc ((nc+2) * sizeof (float)); t0 = (float **) malloc ((nr+2) * sizeof (float *)); for (i = 0; i < nr+2; i++) t0[i] = (float *) malloc ((nc+2) * sizeof (float)); g0 = (float **) malloc ((nr+2) * sizeof (float *)); for (i = 0; i < nr+2; i++) g0[i] = (float *) malloc ((nc+2) * sizeof (float)); int1v = (float **) malloc ((nr+2) * sizeof (float *)); for (i = 0; i < nr+2; i++) int1v[i] = (float *) malloc ((nc+2) * sizeof (float)); int2v = (float **) malloc ((nr+2) * sizeof (float *)); for (i = 0; i < nr+2; i++) int2v[i] = (float *) malloc ((nc+2) * sizeof (float)); int3v = (float **) malloc ((nr+2) * sizeof (float *)); for (i = 0; i < nr+2; i++) int3v[i] = (float *) malloc ((nc+2) * sizeof (float)); leader = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) leader[i] = (int *) malloc ((nc+2) * sizeof (int)); livebrood = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) livebrood[i] = (int *) malloc ((nc+2) * sizeof (int)); newad = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) newad[i] = (int *) malloc ((nc+2) * sizeof (int)); oldad = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) oldad[i] = (int *) malloc ((nc+2) * sizeof (int)); newovi = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) newovi[i] = (int *) malloc ((nc+2) * sizeof (int)); oldovi = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) oldovi[i] = (int *) malloc ((nc+2) * sizeof (int)); newjuv = (int ***) malloc ((nr+2) * sizeof (int **)); for (i = 0; i < nr+2; i++) { newjuv[i] = (int **) malloc ((nc+2) * sizeof (int *)); for (j = 0; j < nc+2; j++) newjuv[i][j] = (int *) malloc (200 * sizeof (int)); } oldjuv = (int ***) malloc ((nr+2) * sizeof (int **)); for (i = 0; i < nr+2; i++) { oldjuv[i] = (int **) malloc ((nc+2) * sizeof (int *)); for (j = 0; j < nc+2; j++) oldjuv[i][j] = (int *) malloc (200 * sizeof (int)); } maxjuv = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) maxjuv[i] = (int *) malloc ((nc+2) * sizeof (int)); emerge = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) emerge[i] = (int *) malloc ((nc+2) * sizeof (int)); newdeg = (float ***) malloc ((nr+2) * sizeof (float **)); for (i = 0; i < nr+2; i++) { newdeg[i] = (float **) malloc ((nc+2) * sizeof (float *)); for (j = 0; j < nc+2; j++) newdeg[i][j] = (float *) malloc (200 * sizeof (float)); } olddeg = (float ***) malloc ((nr+2) * sizeof (float **)); for (i = 0; i < nr+2; i++) { olddeg[i] = (float **) malloc ((nc+2) * sizeof (float *)); for (j = 0; j < nc+2; j++) olddeg[i][j] = (float *) malloc (200 * sizeof (float)); } newh = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) newh[i] = (int *) malloc ((nc+2) * sizeof (int)); oldh = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) oldh[i] = (int *) malloc ((nc+2) * sizeof (int)); neweh = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) neweh[i] = (int *) malloc ((nc+2) * sizeof (int)); oldeh = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) oldeh[i] = (int *) malloc ((nc+2) * sizeof (int)); newvol = (float **) malloc ((nr+2) * sizeof (float *)); for (i = 0; i < nr+2; i++) newvol[i] = (float *) malloc ((nc+2) * sizeof (float)); oldvol = (float **) malloc ((nr+2) * sizeof (float *)); for (i = 0; i < nr+2; i++) oldvol[i] = (float *) malloc ((nc+2) * sizeof (float)); create = (float **) malloc ((nr+2) * sizeof (float *)); for (i = 0; i < nr+2; i++) create[i] = (float *) malloc ((nc+2) * sizeof (float)); repel = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) repel[i] = (int *) malloc ((nc+2) * sizeof (int)); postlarv2 = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) postlarv2[i] = (int *) malloc ((nc+2) * sizeof (int)); numegg = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) numegg[i] = (int *) malloc ((nc+2) * sizeof (int)); totnewegg = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) totnewegg[i] = (int *) malloc ((nc+2) * sizeof (int)); totemerge = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) totemerge[i] = (int *) malloc ((nc+2) * sizeof (int)); numlarv = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) numlarv[i] = (int *) malloc ((nc+2) * sizeof (int)); numpupa = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) numpupa[i] = (int *) malloc ((nc+2) * sizeof (int)); youngad = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) youngad[i] = (int *) malloc ((nc+2) * sizeof (int)); totjuv = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) totjuv[i] = (int *) malloc ((nc+2) * sizeof (int)); randif = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) randif[i] = (int *) malloc ((nc+2) * sizeof (int)); damage = (float **) malloc ((nr+2) * sizeof (float *)); for (i = 0; i < nr+2; i++) damage[i] = (float *) malloc ((nc+2) * sizeof (float)); wmerch = (float **) malloc ((nr+2) * sizeof (float *)); for (i = 0; i < nr+2; i++) wmerch[i] = (float *) malloc ((nc+2) * sizeof (float)); wee_tree = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) wee_tree[i] = (int *) malloc ((nc+2) * sizeof (int)); adjacent = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) adjacent[i] = (int *) malloc ((nc+2) * sizeof (int)); treetype = (int **) malloc ((nr+2) * sizeof (int *)); for (i = 0; i < nr+2; i++) treetype[i] = (int *) malloc ((nc+2) * sizeof (int)); blocktype = (int **) malloc ((nrblk+2) * sizeof (int *)); for (i = 0; i < nrblk+2; i++) blocktype[i] = (int *) malloc ((ncblk+2) * sizeof (int)); clonelist = (int *) malloc ((nrblk*ncblk+1) * sizeof (int *)); clntree = (int *) malloc ((numcln+1) * sizeof (int *)); startclntree = (int *) malloc ((numcln+1) * sizeof (int *)); startw = (float *) malloc ((numcln+1) * sizeof (float *)); /* Call a random seed so the reps will all be different */ /* Note - We don't need this now, as the initial seed values are * already randomized */ InitRandom (firstseed); /* Parameter values for trees chosen for rapid growth */ g0mean = 1.0; g0sd = 0.1; g0var = 0.01; rho = 0.5; r0sd = sqrt (R0VAR); g0r0var = g0var * (1.0 - rho*rho); hsq = 0.5; /* Heritability */ if ( hsq > 0.0 ) sitevar = g0var*(1.0 - hsq)/hsq; /* Variance component of g0 that is * due to site variation. For now, * heritibility is chosen to be 0.5, * so sitevar g0var */ /* When rancln = 0, clone parameters are specified, not random. `High' * values for r0, a0, and t0 are 2 std dev above the mean, and `low' * values are 2 std dev below the mean. Values of g0 are correlated with r0. * Here, high values are 2 std above the conditional mean (conditional on * r0) and low values are std dev below this mean. */ if (rancln == 0) { a0cl[1] = 8.0; r0cl[1] = 2.0; t0cl[1] = 2.0; g0cl[1] = 0.8; a0cl[2] = 2.0; r0cl[2] = 8.0; t0cl[2] = 8.0; g0cl[2] = 1.2; a0cl[3] = 8.0; r0cl[3] = 8.0; t0cl[3] = 2.0; g0cl[3] = 1.2; a0cl[4] = 2.0; r0cl[4] = 2.0; t0cl[4] = 8.0; g0cl[4] = 0.8; a0cl[5] = 2.0; r0cl[5] = 8.0; t0cl[5] = 2.0; g0cl[5] = 1.2; a0cl[6] = 5.0; r0cl[6] = 5.0; t0cl[6] = 5.0; g0cl[6] = 1.0; a0cl[7] = 2.0; r0cl[7] = 2.0; t0cl[7] = 2.0; g0cl[7] = 0.8; a0cl[8] = 8.0; r0cl[8] = 2.0; t0cl[8] = 8.0; g0cl[8] = 0.8; a0cl[9] = 8.0; r0cl[9] = 8.0; t0cl[9] = 8.0; g0cl[9] = 1.2; a0cl[10] = 8.0; r0cl[10] = 2.0; t0cl[10] = 2.0; g0cl[10] = 1.1; a0cl[11] = 2.0; r0cl[11] = 8.0; t0cl[11] = 8.0; g0cl[11] = 0.9; a0cl[12] = 8.0; r0cl[12] = 8.0; t0cl[12] = 2.0; g0cl[12] = 0.9; a0cl[13] = 2.0; r0cl[13] = 2.0; t0cl[13] = 8.0; g0cl[13] = 1.1; a0cl[14] = 2.0; r0cl[14] = 8.0; t0cl[14] = 2.0; g0cl[14] = 0.9; a0cl[15] = 5.0; r0cl[15] = 5.0; t0cl[15] = 5.0; g0cl[15] = 1.0; a0cl[16] = 2.0; r0cl[16] = 2.0; t0cl[16] = 2.0; g0cl[16] = 1.1; a0cl[17] = 8.0; r0cl[17] = 2.0; t0cl[17] = 8.0; g0cl[17] = 1.1; a0cl[18] = 8.0; r0cl[18] = 8.0; t0cl[18] = 8.0; g0cl[18] = 0.9 ; /* Site variation will be added when ramets are planted */ } /* End if rancln = 0 */ if (rancln == 1) /* Clone parameters are generated randomly */ { /* * Generate the clones to be used (total number is numcln) * in planting the stand. Clones have minimum toxicity, * attractiveness, and resin flow rates (t0, a0, r0, respectively) * of 0.1, but there are no upper bounds on these quantities. * Tree vigour g0 is correlated with r0. */ for (n=1; n<=numcln; n++) { t0cl[n] = NormalRandom(T0MEAN,T0VAR); a0cl[n] = NormalRandom(A0MEAN,A0VAR); r0cl[n] = NormalRandom(R0MEAN,R0VAR); if (a0cl[n]<0.1) a0cl[n]=0.1; if (r0cl[n]<0.1) r0cl[n]=0.1; if (t0cl[n]<0.1) t0cl[n]=0.1; g0r0mean = g0mean + (rho * g0sd / r0sd) * (r0cl[n] - R0MEAN); g0cl[n] = NormalRandom (g0r0mean, g0r0var); if (g0cl[n]<0.1) g0cl[n]=0.1; } /* End n loop */ /* Site variation will be added when ramets are planted */ } /* End if rancln = 1 */ /* 'PLANT' THE STAND AND LAY OUT THE INITIAL WEEVIL DISTRIBUTION */ /* PLANTING THE STAND -- First, zero out the border */ for (i=0; i<=nr+1; i++) { a0[i][0]=0.0; /* NB -- Need a0 = 0 on the border */ r0[i][0]=0.0; /* In the code, r0, g0 and t0 actually don't */ t0[i][0]=0.0; /* matter on the border */ g0[i][0]=0.0; a0[i][nc+1]=0.0; r0[i][nc+1]=0.0; t0[i][nc+1]=0.0; g0[i][nc+1]=0.0; wee_tree[i][0] = 0; /* There are trees only in the interior */ wee_tree[i][nc+1] = 0; /* not on the border */ } /* End i loop */ #ifdef DEBUG printf ("weevils.c: weevil_plant(): done i loop\n"); #endif for (j=1; j<=nc; j++) { a0[0][j]=0.0; r0[0][j]=0.0; t0[0][j]=0.0; g0[0][j]=0.0; a0[nr+1][j]=0.0; r0[nr+1][j]=0.0; t0[nr+1][j]=0.0; g0[nr+1][j]=0.0; wee_tree[0][j] = 0; /* There are trees only in the interior */ wee_tree[nr+1][j] = 0; /* not on the border */ for (i=1; i<=nr; i++) wee_tree[i][j] = 1; /* There ARE trees in the interior */ } /* End j loop */ #ifdef DEBUG printf ("weevils.c: weevil_plant(): done j loop\n"); #endif /* NOTE -- The user can specify a desired planting pattern by inserting his or her own code here. Otherwise, one of the following alternatives can be chosen. */ /* HERE ARE THE ALREADY CODED PATTERNS */ /* Pattern = 0 means each tree is a different 'clone' */ if ( pattern == 0 ) { for (i=1 ; i<=nr ; i++) { for (j=1 ; j<=nc ; j++) { a0[i][j] = NormalRandom(A0MEAN,A0VAR) ; r0[i][j] = NormalRandom(R0MEAN,R0VAR) ; t0[i][j] = NormalRandom(T0MEAN,T0VAR) ; if (a0[i][j] < 0.1) a0[i][j] = 0.1; if (r0[i][j] < 0.1) r0[i][j] = 0.1; if (t0[i][j] < 0.1) t0[i][j] = 0.1; g0r0mean = g0mean + (rho * g0sd / r0sd) * (r0[i][j] - R0MEAN); /* Compute a value of g0, conditional on r0 */ g0[i][j] = NormalRandom (g0r0mean, g0r0var); if (g0[i][j] < 0.1) g0[i][j] = 0.1; /* Add site variation */ g0[i][j] = g0[i][j] + NormalRandom (0.0, sitevar); if (g0[i][j] < 0.1) g0[i][j] = 0.1; itn = (i-1)*nc+(j-1)+1; tree[itn]->data->ht_vig = g0[i][j] ; if (tree[itn]->data->ht_vig < 0.1) tree[itn]->data->ht_vig = 0.1; } /* End for j */ } /* End for i */ } /* End if pattern = 0 */ /* When pattern = 11, trees are planted in a random mix of already selected clones */ if (pattern == 11) { #ifdef DEBUG printf ("weevils.c: weevil_plant(): pattern=11\n"); #endif for (i=1; i<=nr; i++) { for (j=1; j<=nc; j++) { /* urand is a uniform random number distributed U(0,1) */ urand=UniformRandom() ; randcl=ceil(urand*numcln) ; treetype[i][j] = randcl ; a0[i][j]=a0cl[randcl] ; r0[i][j]=r0cl[randcl] ; t0[i][j]=t0cl[randcl] ; g0[i][j]=g0cl[randcl] ; /* Add site variation */ g0[i][j] = g0[i][j] + NormalRandom (0.0, sitevar); if (g0[i][j] < 0.1) g0[i][j] = 0.1; itn = (i-1)*nc+(j-1)+1; tree[itn]->data->ht_vig = g0[i][j]; if (tree[itn]->data->ht_vig < 0.1) tree[itn]->data->ht_vig = 0.1; } /* End j loop */ } /* End i loop */ #ifdef DEBUG printf ("weevils.c: weevil_plant(): done pattern 11\n"); #endif } /* End pattern = 11 */ if ( pattern == 12 || pattern == 183 || pattern == 303 ) { /* Plant already selected clones in randomly chosen blocks; * each clone is used bpc times. Pattern 12 can be used to * designate any such pattern except 183 or 303 */ /* clonelist[k] contains the number of the clone in block ij, * where k = (i-1)*ncblk + j; */ for (i=1; i<=numcln; i++) { for (j=1; j<=bpc; j++) { k = (i-1)*bpc + j; clonelist[k] = i; } /* End j */ } /* End i */ for (k=1; k<=numcln*bpc; k++) { urand = UniformRandom() ; randcl = ceil(urand*numcln*bpc) ; num2 = clonelist[k]; clonelist[k] = clonelist[randcl]; clonelist[randcl] = num2; } /* End k */ for (k=1; k<=nrblk; k++) { for (m=1; m<=ncblk; m++) { num1 = (k-1)*rpb+1; num2 = k*rpb; num3 = (m-1)*cpb+1; num4 = m*cpb; n = (k-1)*ncblk + (m-1)+1; blocktype[k][m] = clonelist[n] ; for (i=num1; i<=num2; i++) { for (j=num3; j<=num4; j++) { a0[i][j] = a0cl[clonelist[n]] ; r0[i][j] = r0cl[clonelist[n]] ; t0[i][j] = t0cl[clonelist[n]] ; g0[i][j] = g0cl[clonelist[n]] ; treetype[i][j] = clonelist[n]; /* Add site variation */ g0[i][j] = g0[i][j] + NormalRandom (0.0, sitevar); itn = (i-1)*nc+(j-1)+1; tree[itn]->data->ht_vig = g0[i][j]; if (tree[itn]->data->ht_vig < 0.1) tree[itn]->data->ht_vig = 0.1; } /* End j loop */ } /* End i loop */ } /* End m loop */ } /* End k loop */ } /* End if pattern = 12 */ /* SOME FIXED PATTERNS USING 2, 6, 18 or 30 CLONES */ if ( pattern == 21 ) { clonenumber[1] = 1; clonenumber[2] = 2; } /* End if pattern = 21 */ if (pattern == 22 ) { clonenumber[1] = 1; clonenumber[2] = 2; clonenumber[3] = 2; clonenumber[4] = 1; } /* End if pattern = 22 */ if (pattern == 23) { clonenumber[1] = 1; clonenumber[2] = 2; clonenumber[3] = 1; clonenumber[4] = 2; clonenumber[5] = 2; clonenumber[6] = 1; clonenumber[7] = 2; clonenumber[8] = 1; clonenumber[9] = 1; clonenumber[10] = 2; clonenumber[11] = 1; clonenumber[12] = 2; clonenumber[13] = 2; clonenumber[14] = 1; clonenumber[15] = 2; clonenumber[16] = 1; } /* End if pattern = 23 */ if ( pattern == 61 ) { clonenumber[1] = 1; clonenumber[2] = 2; clonenumber[3] = 3; clonenumber[4] = 4; clonenumber[5] = 5; clonenumber[6] = 6; } /* End if pattern = 61 */ if (pattern == 62) { clonenumber[1] = 1; clonenumber[2] = 6; clonenumber[3] = 5; clonenumber[4] = 4; clonenumber[5] = 3; clonenumber[6] = 2; clonenumber[7] = 2; clonenumber[8] = 4; clonenumber[9] = 6; clonenumber[10] = 3; clonenumber[11] = 1; clonenumber[12] = 5; clonenumber[13] = 5; clonenumber[14] = 1; clonenumber[15] = 4; clonenumber[16] = 2; clonenumber[17] = 6; clonenumber[18] = 3; clonenumber[19] = 4; clonenumber[20] = 3; clonenumber[21] = 2; clonenumber[22] = 1; clonenumber[23] = 5; clonenumber[24] = 6; clonenumber[25] = 6; clonenumber[26] = 2; clonenumber[27] = 3; clonenumber[28] = 5; clonenumber[29] = 4; clonenumber[30] = 1; clonenumber[31] = 3; clonenumber[32] = 5; clonenumber[33] = 1; clonenumber[34] = 6; clonenumber[35] = 2; clonenumber[36] = 4; } /* End if pattern = 62 */ if (pattern == 181 || pattern == 301 ) { for (k=1; k<=numcln; k++) clonenumber[k] = k; } /* End if pattern = 181 or pattern = 301 */ /* PLANT TREES IN BLOCKS USING THE ABOVE SPECIFIED BLOCK PATTERNS */ if ( pattern != 0 && pattern != 11 && pattern != 12 && pattern != 183 && pattern != 303) { for (k=1; k<=nrblk; k++) { for (m=1; m<=ncblk; m++) { num1 = (k-1)*rpb+1; num2 = k*rpb; num3 = (m-1)*cpb+1; num4 = m*cpb; n = (k-1)*ncblk + (m-1)+1; blocktype[k][m] = clonenumber[n] ; for (i=num1; i<=num2; i++) { for (j=num3; j<=num4; j++) { a0[i][j] = a0cl[clonenumber[n]] ; r0[i][j] = r0cl[clonenumber[n]] ; t0[i][j] = t0cl[clonenumber[n]] ; g0[i][j] = g0cl[clonenumber[n]] ; treetype[i][j] = clonenumber[n] ; /* Add site variation */ g0[i][j] = g0[i][j] + NormalRandom (0.0, sitevar); itn = (i-1)*nc+(j-1)+1; tree[itn]->data->ht_vig = g0[i][j]; if (tree[itn]->data->ht_vig < 0.1) tree[itn]->data->ht_vig = 0.1; } /* End j loop */ } /* End i loop */ } /* End m loop */ } /* End k loop */ } /* End if pattern != 0, etc. */ /* THIS ENDS PLANTING OF THE STAND OF TREES */ /* TWO CONSTANTS RELATED TO DIFFUSION OF VOLATILES AND WEEVILS */ vstay = 1.0 - exp(-1.7*1.7/(4.0*RVDIFFUS)); /* * The particular form of vstay reflects the fact that the pattern of * diffused particles originating at a point is that of a 2-D normal * density */ wstay = 1.0 - exp(-1.7*1.7/(4.0*RADDIFFUS)); /* * The one-day time interval is short enough that we * can assume that if weevils diffuse, they move no more than * one tree (grid square) away. */ tol = 0.01*nr*nc ; /* Relates to convergence of the recursion for volatiles */ /* GET THE SIMULATION STARTED - ZERO OUT SOME MATRICES */ for (i=0; i<=nr+1; i++) { for (j=0; j<=nc+1; j++) { newvol[i][j] = 1.0; /* Just to get started */ livebrood[i][j] = 0; newad[i][j] = 0; crown[i][j] = 0.0; length[i][j] = 0.0; height[i][j] = 0.0; diam[i][j] = 0.0; create[i][j] = 0.0; randif[i][j] = 0; } /* End for j */ } /* End for i */ /* FREE ARRAYS NO LONGER NEEDED */ #ifdef DEBUG printf ("weevils.c: weevil_plant(): free memory\n"); #endif free (clonelist); for (i = 0; i < nrblk+2; i++) free (blocktype[i]); free (blocktype); /* reset random number seed */ /* Note - We don't need this now, as the initial seed values are * already randomized */ InitRandom (secondseed); /* Count the stems planted of each clone type */ if (pattern != 0) { for (m=1; m<=numcln; m++) { startclntree[m] = 0; for (i=1; i<=nr; i++) { for (j=1; j<=nc; j++) { if (treetype[i][j] == m) startclntree[m] = startclntree[m] + wee_tree[i][j]; } /* end j loop */ } /* end i loop */ } /* End m loop */ } /* End if pattern != 0 */ #ifdef DEBUG printf ("weevils.c: weevil_plant(): end\n"); #endif } /* end of weevil_plant() */ /* Stand planting ends; simulation begins */ /* The loop on nyear starts here */ void weevil_year (FILE *errfile, int nyear, struct weevil_parameters *params) { /* Vectors related to printing of output files */ float sumdcl[51]; float percentcl[51]; float summerchcl[51]; float avhgtcl[51]; int sumattkcl[51]; /* total number of attacked trees in each clone */ float percentattkcl[51]; /* percent attacked trees in each clone */ /* NONOVIPOSITORS */ float wintersurv ; /* adult over-winter survival rate */ /* PARAMETER IN THE EQUATION FOR EGG HOLES AND EGGS */ float re1; /* rate at which individual ovipositing weevils * deposit eggs. The value of this parameter depends on * temperature and is computed in the code */ /* WIND */ float wind[2]; /* wind = 2-D wind vector: wind(0) is east-west, * positive direction is east; wind(1) is n-s * pos. direction is south */ float unit; /* Temporary variable re wind vector */ float wspeed; /* wind speed, randomly generated each day */ /* TEMPERATURE */ float totdeg ; /* total number of degree days in a year */ float theta; float temp; /* today's 'actual' average temperature - a normal * random variable, computed */ float tmin; /* today's minimum temperature ( = temp - tdev) */ float tdmean; /* mean of tdev - about 5 degrees C in BC locations; this varies over the year */ float tdev; /* magnitude of temperature deviation around temp - a gamma r.v. different mean and variance in winter than in summer */ float xdeg; /* number of degree-days accumulated 'today' - randomly * generated using the following quantities */ float tmax; /* today's maximum temperature ( = temp + tdev) */ float tavg; /* long-term average temp for 'today' - from weather records, * but approximated by a sine curve */ /* VARIABLES USED IN I,J LOOPS FOR ODORS AND VOLATILES */ float err; float cnsum; /* VARIABLES THAT INDICATE THE STATUS OF TREES AND THEIR WEEVIL POPULATIONS */ /* Variables and vectors relating to average numbers of weevils per tree */ float avnonovi; /* Average number of adult nonovipositors per live tree */ float avovi; /* Average number of ovipositors per live tree */ float avegg; /* Average number of weevil eggs per live tree */ float avlarv; /* Average number of larvae per live tree */ float avpupa; /* Average number of pupae per live tree */ float avyoungad; /* Average number of young (immature) adults per live tree */ float avnewegg; /* Average total number of eggs laid per attacked tree, * this YEAR. Other numbers above are averages as of TODAY */ float avemerge; /* Average total number of emerged brood adults per * attacked tree, this year. */ float avlivebrood; /* Average total number of surviving brood adults per * attacked tree, this year. */ float attract[3][3]; float percent; /* overall percentage of leaders that survive 'this season' */ float weestart; /* overall average no. weevils per tree at start of a season */ int numtrees; /* total number of live trees in the plantation */ float death; /* Juvenile death rate */ float dstuck; /* dstuck = probability of getting stuck and dying while * trying to emerge - emergence is easier, and thus this * probability is smaller, when there are other young * adults that can help pave the way */ int diff; /* temporary variable re juveniles */ int nonovi; /* temporary value for number of non-ovipositing adults */ float term; float xdays; int xnum; float xlength; /* short-hand access to length[][] */ float xa0; /* short-hand access to a0[][] */ float xr0; /* short-hand access to r0[][] */ int xad; /* short-hand access to oldad[][] */ int xh; /* short-hand access to oldh[][] */ float xt0; /* short-hand access to t0[][] */ int xovi; /* short-hand access to oldovi[][] */ int xeh; /* short-hand access to oldeh[][] */ float xcrown; /* short-hand access to crown[][] */ float xdiam; /* short-hand access to diam[][] */ int itn; /* TASS tree number */ float sumd; int sumattk; /* total number of attacked trees in the stand */ float percentattk; /* percent attacked trees in the stand */ float avhgt; float yyy; int nyy, nzz; float sum; int nday; /* day of year (1-365) */ float factor; int num1, i, j, k, n, m; int numwee; /* Initial number of weevils on each tree when pattern = 1 */ int wee2; /* Defined re weedist below */ /* used for computing TASS merch volume */ struct merch_record mrec; #ifdef DEBUG printf ("weevils.c: weevil_year(): begin year %d\n", nyear); #endif if (nyear == 1) weevil_plant (errfile, params); if (nyear < 9) { #ifdef DEBUG printf ("weevils.c: weevil_year(): end\n"); #endif return; } /* End if nyear < 9 */ if (nyear == 9) { /* Specify the initial weevil distribution in the plantation */ if (weedist == 1) { if (noweevil == 1) numwee = 0; else numwee = 3; for (i=1; i<=nr; i++) { for (j=1; j<=nc; j++) { itn = (i-1)*nc+(j-1)+1; if (tree_alive (itn) == TRUE) newad[i][j] = numwee ; } /* End for j */ } /* End for i */ } /* End if weedist = 1; i.e., end INITDIST1 */ if (weedist == 2) { #ifdef DEBUG printf ("weevils.c: weevil_plant(): weedist=2\n"); #endif if (noweevil == 1) wee2 = 0; else wee2 = 20; for (i=1; i<=5; i++) { for (j=1; j<=5; j++) { itn = (i-1)*nc+(j-1)+1; if (tree_alive (itn) == TRUE) newad[i][j] = wee2 ; } /* End for j */ } /* End for i */ } /* End if weedist = 2; i.e., end INITDIST2 */ } /* End if nyear = 9 */ /* INITIAL CONDITIONS FOR THE SUMMER SIMULATION * Start the simulation on March 1, just as the temperature is * high enough so that degree-days begin to accumulate regularly. * When winter temps are simulated and used, this will be deleted. */ tavg = 5.75; /* Average temp in Vancouver on March 1 */ temp = tavg; /* Just to get the simulation started */ totdeg = 0.0; /* total # degree-days this year */ sumd = 0.0; /* Tree related values - Remain zero until day 365 */ sumattk = 0; percent = 0.0; summerch = 0.0; avhgt = 0.0; /* Zero out some matrices */ for (m=1; m<=51; m++) { sumdcl[m] = 0.0; sumattkcl[m] = 0; percentcl[m] = 0.0; summerchcl[m] = 0.0; avhgtcl[m] = 0.0; /* End initialization of tree values */ } /* End loop on m */ /* Initialize weevilvalues */ for (i=1; i<=nr; i++) { for (j=1; j<=nc; j++) { livebrood[i][j] = 0; newovi[i][j] = 0; newh[i][j] = 0; neweh[i][j] = 0; maxjuv[i][j] = 0; emerge[i][j] = 0; damage[i][j] = 0.0; totnewegg[i][j] = 0; totemerge[i][j] = 0; for (n = 0; n<=198; n++) { /* initially zero; values will be added * later as the simulation proceeds */ newdeg[i][j][n] = 0.0 ; newjuv[i][j][n] = 0 ; } /* End for n */ /* In the next section, variables come from TASS * wee_tree[i][j] = 1 if tree ij is alive, and = 0 if it is dead * CROWN is measured in cubic meters * LENGTH is measured in centimeters * DIAMETER is measured in centimeters */ #ifdef DEBUG printf ("weevils.c: weevil_year(): before TASS stuff\n"); #endif itn = (i-1)*nc+(j-1)+1; #ifdef DEBUG printf ("weevils.c: weevil_year(): itn=%d, num_trees=%d\n", itn, num_trees); #endif if (tree_alive (itn) == TRUE) { crown[i][j] = tree[itn]->data->fv; /* Leader length */ length[i][j] = 100*(tree[itn]->data->ht_inc[0]); if (length[i][j] < 0.1) length[i][j] = 0.1 ; if (crown[i][j] < 0.00001) crown[i][j] = 0.00001 ; diam[i][j] = 0.05*length[i][j] ; /* Mark tree's leader as not attacked */ tree[itn]->data->itemp = 0; leader[i][j] = 1; } else { crown[i][j] = 0.0; height[i][j] = 0.0; length[i][j] = 0.0; diam[i][j] = 0.0; a0[i][j] = 0.0; create[i][j] = 0.0; newvol[i][j] = 0.0; newad[i][j] = 0; wee_tree[i][j] = 0; /* TASS says tree is dead */ leader[i][j] = 0; } /* End else tree is dead */ #ifdef DEBUG printf ("weevils.c: weevil_year(): done TASS\n"); #endif create[i][j] = RVCREATE * crown[i][j] ; } /* End j */ } /* End i */ /* numtrees is the total number of live trees in the stand */ /* Compute the overall average number of weevils per live tree, and the * average number per tree for each clone, at the start of the season */ numtrees = 0; weestart = 0.0 ; for (i=1; i<=nr; i++) { for (j=1; j<=nc; j++) { numtrees = numtrees + wee_tree[i][j]; weestart = weestart + newad[i][j] ; } /* end j loop */ } /* end i loop */ if ( numtrees > 0) weestart = weestart/numtrees ; else weestart = 0.0 ; if (pattern != 0) { for (m=1; m<=numcln; m++) { startw[m] = 0.0; clntree[m] = 0; for (i=1; i<=nr; i++) { for (j=1; j<=nc; j++) { if (treetype[i][j] == m) { startw[m] = startw[m] + newad[i][j] ; clntree[m] = clntree[m] + wee_tree[i][j]; } /* wee_tree and newad are zero if a tree is dead */ } /* end j loop */ } /* end i loop */ if (clntree[m] != 0) startw[m] = startw[m]/clntree[m] ; else startw[m] = 0.0 ; } /* End m loop */ } /* End if pattern != 0 */ /* Compute matrix 'adjacent' : * adjacent[i][j] is the number of live trees adjacent to tree ij */ #ifdef DEBUG printf ("weevils.c: weevil_year(): before computing matrix adjacent\n"); #endif for (i=1; i<=nr; i++) { for (j=1; j<=nc; j++) { adjacent[i][j] = 0 ; for ( k=0 ; k<=2 ; k++ ) { for ( m=0 ; m<=2 ; m++ ) adjacent[i][j] = adjacent[i][j] + wee_tree[i-1+k][j-1+m] ; /* This computation works because entries in matrix 'wee_tree' * are zero on the border */ } /* End for k */ adjacent[i][j] = adjacent[i][j] - wee_tree[i][j] ; /* Tree i,j is not adjacent to itself */ } /* end j loop */ } /* end i loop */ /* set up for computing TASS merch volume per tree */ mrec.top_dib = 10; mrec.min_dbh = 12.5; mrec.bot_dib = 0; mrec.bot_dob = 0; mrec.stump_height = 0.30; /* WITHIN EACH YEAR SIMULATE DAILY WEEVIL POPULATION CHANGES */ /* If noweevil = 1, there are never any weevils, and only TASS (tree) * computations produce effects on the stand */ /* The loop on nday starts here */ for (nday=60; nday<=345; nday++) { #ifdef DEBUG printf ("weevils.c: weevil_year(): nday=%d\n", nday); #endif /* * GENERATE TODAY'S WIND VELOCITY VECTOR * * Wind speed ranges from 0 to 2 (gridlines/sec, about 0 to 14mph), * with a 100*WEW % chance of blowing from the west and a 100*WNS % * chance of blowing from the south", wind[0] blows east-west, * while wind[1] blows north-south * * First, establish the DIRECTIONS - Speed will be determined later */ urand = UniformRandom() ; wind[0] = WEW - urand ; if (wind[0] > 0.0) wind[0] = wind[0] / WEW ; else wind[0] = wind[0] /(1.0 - WEW); urand = UniformRandom() ; wind[1] = WNS - urand ; if (wind[1] > 0.0) wind[1] = wind[1] /(1.0 - WNS) ; else wind[1] = wind[1] / WNS ; /* * The effect is that wind[0] or wind[1] can be in either direction, * and can have any length between 0 and 1 */ xxx = wind[1] * wind[1] + wind[0] * wind[0]; if ( xxx <= 0.0) unit = 0.0 ; else unit = 1.0/(sqrt(xxx)) ; wind[1] = wind[1] * unit ; wind[0] = wind[0] * unit ; /* * At this point, wind[0] and wind[1] constitute a unit vector in the * wind direction. Now impose wind speed. */ urand = UniformRandom() ; wspeed = 22.0*urand ; /* * In this version, wind speed is a uniform random variable; later, I'll * make this more realistic by changing to a gamma distribution, * in which there is no max speed, but speeds are more likely to be * small than large. For now, 22.0 is the maximum wind * speed, in km/hr. Now change to grid lines per second. * With the current 3 meter grid, max wind speed is about * 2 grid lines per second (about 14 mph). */ wspeed = 10.0*wspeed/(36.0*GRID) ; wind[1] = wind[1] * wspeed ; wind[0] = wind[0] * wspeed ; /* COMPUTE TODAY'S TEMPERATURES AND DEGREE-DAYS */ tavg = tmean - TAMPL*cos(2.0*PI*(nday-20.0)/365.0); /* * Today's long-term avg temp, tavg, follows a cosine curve over the * year. A year contains 365 days. Day 1 is Jan 1; Jan 20 is the * coldest day. */ tavg = 0.5*(tavg + temp); /* * The long-term average is replaced by the average of the true * long-term average and yesterday's observed average, temp, thus * correlating successive daily temps. * * temp is then updated to equal today's actual (i.e., simulated) average * temperature, chosen randomly from a normal distribution having * mean = tavg */ temp = NormalRandom(tavg,TVAR); /* tdmean = avg difference between daily average and max temperatures * tdev = simulated difference between these two */ tdmean = 3.5 + 1.5*sin(2.0*PI*(nday-110.0)/365.0); tdev = fabs(NormalRandom(tdmean, TDVAR)); /* tmax is today's simulated maximum temperature */ tmax = temp + tdev ; tmin = temp - tdev ; /* not used right now * tlambda = tdmean/TDVAR ; * tee = tlambda* tdmean ; */ if ( tmax < THRESH ) xdeg = 0.0 ; else if ( tmin > THRESH ) xdeg = temp - THRESH ; else { xxx = (THRESH - temp)/tdev; if ( xxx > 1.0) xxx = 1.0 ; else if ( xxx < -1.0) xxx = -1.0 ; theta = asin(xxx) ; xdeg = (1.0/PI)*( (temp - THRESH)*((PI/2) - theta) + tdev*cos(theta)); } /* Formula for xdeg from Baskerville and Emin (1969) - Has been checked */ if (xdeg > 0.0) totdeg = totdeg + xdeg ; #ifdef DEBUG printf ("weevils.c: weevil_year(): done temperature calcs\n"); #endif /* UPDATE YESTERDAY'S 'OLD' MATRICES; REPLACE BY TODAY'S NEW ONES */ for (i=1; i<=nr; i++) { for (j=1; j<=nc; j++) { /* Assign mortality to this year's brood adults. * For simplicity, we pretend they do not join the new * adults until the end of the season. Hence, they * do not move from tree to tree, oviposit, or otherwise * contribute to the population until the following spring. * For this year, we assume they just fatten up in order * to survive the coming winter. */ livebrood[i][j] = chance_check(livebrood[i][j], 1.0-DADULT); /* Do these next things only if the tree is alive */ if (wee_tree[i][j] == 0) continue; oldad[i][j] = newad[i][j] ; oldovi[i][j] = newovi[i][j] ; oldeh[i][j] = neweh[i][j] ; oldh[i][j] = newh[i][j] ; for (n = 0; n<=198; n++) { olddeg[i][j][n] = newdeg[i][j][n] ; oldjuv[i][j][n] = newjuv[i][j][n] ; } /* End for n */ /* * Use the 'yesterday' (old) matrices to compute new weevil * populations 'today' */ /* First, shorten the notation */ xad = oldad[i][j]; xovi = oldovi[i][j]; xeh = oldeh[i][j]; xh = oldh[i][j]; xr0 = r0[i][j]; xa0 = a0[i][j]; xt0 = t0[i][j]; xcrown = crown[i][j]; xlength = length[i][j]; xdiam = diam[i][j]; /* * Now for the computations (the daily update). * Start with the things that do not move from tree to tree. */ /* THE RECURSION FOR NUMBER OF EGG HOLES */ if (temp < 10.0 || temp > 35.0 ) re1 = 0.0; else { if (temp > 10.0 && temp < 29.0) re1 = 0.17*exp(0.15*(temp - 10.0)); else re1 = 0.56*(35.0 - temp); } if (re1 < 0.0) re1 = 0.0; re1 = multre1*re1; nyy = PoissonRandom((re1/1.4)*xovi); neweh[i][j] = xeh + nyy; /* * re1 = daily rate at which EGGS are produced by a single ovipositor * there are, on average, 1.4 eggs per egg hole */ /* THE RECURSION FOR TOTAL NUMBER OF HOLES */ xxx = xdiam*xlength ; term = 1.0 - RH5*xh/(PI*xxx) ; if (term<0.01) term = 0.01 ; nzz = PoissonRandom(RH1*(xad - xovi)*(exp(-RH3*xr0-RH4*xt0))*sqrt(term)); newh[i][j] = xh + nyy + nzz; /* For ovipositors, egg holes are also feeding holes */ /* THE RECURSION FOR JUVENILES */ { /* * Juvenile action takes place on all days #60-345 having high * enough temps */ /* First, check whether any degree-days were added today */ if (xdeg == 0.0) { /* * Here, number of deg-days = 0 today, and survival is all that * needs to be updated. * * These next computations will be useful for the Weevil Dynamics ms */ postlarv2[i][j] = 0; /* * postlarv2 = yesterday's number of juveniles in * categories larvae 3 and older; i.e., those having * at least 283 degree-days */ for (n = 0; n<=198; n++) { xnum = oldjuv[i][j][n] ; xdays = olddeg[i][j][n]; if (xdays > 283.0 ) postlarv2[i][j] = postlarv2[i][j] + xnum ; } /* end for n */ /* * now have the total # of postlarv2 juveniles from yesterday; * also have the total numbers of all juveniles */ for (n = 0; n<=198; n++) { xnum = oldjuv[i][j][n] ; xdays = olddeg[i][j][n]; /* death = mortality rate */ if (xdays < 142.0) death = 0.01 ; else if (xdays < 283.0 ) death = 0.012 ; else if (xdays < 568.0 ) { /* * Impose mortality due to crowding; the next lines * allow 7 mm of leader circumference per (Check whether we * juvenile before crowding begins to affect want to allow * the mortality rate (See, e.g., Silver, 1968) more larvae per cm) */ xxx = postlarv2[i][j]/(14.0*PI*diam[i][j]); if (xxx > 1.0) death = 0.02*(1.0 + fabs(1.0 - xxx)); else death = 0.02; } else death = 0.02; if (death > 1.0) death = 1.0; /* Now have the death rate appropriate to this juvenile group */ newjuv[i][j][n] = chance_check(xnum, 1.0-death); } /* end for n */ } /* end if deg-days = 0 */ else { /* * Now, number of degree-days is positive. In this case, * the only immature beetles having no degree-days are * new eggs */ emerge[i][j] = 0; newjuv[i][j][0] = PoissonRandom(re1*xovi); totnewegg[i][j] = totnewegg[i][j] + newjuv[i][j][0]; /* postlarv2 = number of juveniles having at least 283 * degree-days; i.e., those in larval instars 3 and 4, * pupae and young adults * * These next computations will be useful for the Weevil Dynamics ms */ postlarv2[i][j] = 0; youngad[i][j] = 0; /* * postlarv2 = yesterday's number of juveniles in categories * larvae 3 and older; i.e., those having at least 283 degree-days * youngad = number of juveniles with at least 737 degree-days */ for (n = 0; n<=198; n++) { xnum = oldjuv[i][j][n] ; xdays = olddeg[i][j][n]; if (xdays > 736.0) youngad[i][j] = youngad[i][j] + xnum; if (xdays > 283.0) postlarv2[i][j] = postlarv2[i][j] + xnum ; } /* end for n */ /* * Now have total numbers of eggs, larvae, pupae, young adults * and postlarv2, as well as the total number of all juveniles, * from yesterday */ diff = postlarv2[i][j] - youngad[i][j]; if (diff > maxjuv[i][j]) maxjuv[i][j] = diff; /* * maxjuv = the largest number in larval stages 3 and 4 in the * feeding ring at any one time this year */ for (n = 0; n<=198; n++) { xnum = oldjuv[i][j][n] ; xdays = olddeg[i][j][n]; /* death = mortality rate */ if (xdays < 142.0) death = 0.01 ; else if (xdays < 283.0 ) death = 0.012 ; else if (xdays < 568.0 ) { /* * Impose mortality due to crowding; the next lines allow 7 mm * of leader circumference per juvenile before crowding begins * to affect the mortality rate (See, e.g., Silver, 1968) */ xxx = postlarv2[i][j]/(14.0*PI*diam[i][j]); if (xxx > 1.0) death = 0.02*(1.0 + fabs(1.0 - xxx)); else death = 0.02; } else death = 0.02; if (death > 1.0) death = 1.0; /* Now have the death rate appropriate to this juvenile group */ newjuv[i][j][n+1] = chance_check(xnum, 1.0-death); newdeg[i][j][n+1] = xdays + xdeg ; /* Compute number that emerge -- For now, make it simple, * but later, make time to full maturity be a * Normal(888,400) random variable */ dstuck = 0.1* exp(-0.1 * (youngad[i][j] - 10.0)) ; if (dstuck > 1.0) dstuck = 1.0 ; if (newdeg[i][j][n+1] > 888.0) { num1 = chance_check(newjuv[i][j][n+1], 1.0 - dstuck); emerge[i][j] = emerge[i][j] + num1; totemerge[i][j] = totemerge[i][j] + num1; livebrood[i][j] = livebrood[i][j] + num1; newdeg[i][j][n+1] = 0.0; newjuv[i][j][n+1] = 0; } /* End IF newdeg > 888 */ } /* end n loop */ } /* end ELSE xdeg > 0 */ } /* end if nday re computation for juveniles */ /* THE RECURSION FOR NUMBER OF OVIPOSITING WEEVILS */ /* Females only -- Assume a 50-50 sex ratio */ if ( (nday <= 59) || (nday >= 346) ) { newovi[i][j] = 0; newad[i][j] = chance_check(xad, 1.0 - DADULT); /* Mortality of adults */ } else { /* Oviposition takes place whenever the temperature is high enough */ nonovi = xad - xovi; nonovi = chance_check(nonovi, 1.0 - DADULT); /* Mortality of nonovipositors */ xovi = chance_check(xovi, 1.0 - DOVI); /* Mortality of ovipositors */ newad[i][j] = xovi + nonovi; xad = newad[i][j]; /* * Next, subtract those that cease oviposition, and add those * that begin; then, make sure the result is not negative */ xxx = xdiam*xlength ; factor = xh/(PI*xxx) ; if (factor < 0.01) factor = 0.01 ; else if (factor > 100.0) factor = 100.0 ; term = 1.0 - RH5*factor ; if (term < 0.01) term = 0.01 ; xxx = ROV1*(1.0+ROV2*xr0+ROV3*xt0)*(1.0-ROV4+ROV4/sqrt(term)); yyy = RNOV1*exp(-RNOV2*xr0 - RNOV3*xt0)*exp(-RNOV4*factor); if(xxx < 0.0) xxx = 0.0; if(xxx > 1.0) xxx = 1.0; if(yyy < 0.0) yyy = 0.0; if(yyy > 1.0) yyy = 1.0; newovi[i][j] = xovi - chance_check(xovi, xxx) + chance_check(nonovi, yyy); /* Adults INCLUDE ovipositors, so no adult adjustment needed; however, adults also include males, so we want to restrict the proportion of ovipositors in the overall population */ } /* end else nday re computation for ovipositors */ } /* End j loop */ } /* End i loop */ /* * NOW FOR VOLATILES AND ADULT WEEVILS, THE THINGS THAT 'MOVE' * THE RECURSION FOR CONCENTRATION OF VOLATILES * * Moving weevils constitute the only category affected by volatile * concentrations. Thus the placement of this section of the code. * Details re volatiles are different from those in all other variable * updates because volatile processes take place on a faster time * scale (seconds) than those of the other variables (days). What * this code does is establish an equilibrium value for volatile * concentration first, then have moving weevils react to that. */ /* THE RECURSION FOR VOLATILES */ err = nr*nc; for ( i=0 ; i<=nr+1 ; i++ ) { for ( j=0 ; j<=nc+1 ; j++ ) oldvol[i][j] = newvol[i][j]; } /* End i loop */ n = 0 ; while (err > tol) /* the tolerance may be changed here */ { /* * DECOMPOSITION -- Multiply preceding volatile matrix by * decomposition rate. * int2v is the matrix of volatiles after decomposition. */ for ( i=0 ; i<=nr+1 ; i++ ) { for ( j=0 ; j<=nc+1 ; j++ ) int2v[i][j] = RVDECOMP*oldvol[i][j] ; } /* End i loop */ /* * That finishes decomposition. Now for DIFFUSION; * int3v is the matrix of volatiles after diffusion. * First, interiors of plantation borders */ for ( i=1; i<=nr ; i++ ) { int3v[i][0] = vstay*int2v[i][0] + ((1-vstay)/8)*(int2v[i-1][0] + int2v[i+1][0] + int2v[i-1][1] + int2v[i][1] + int2v[i+1][1]); int3v[i][nc+1] = vstay*int2v[i][nc+1] + ((1-vstay)/8)*(int2v[i-1][nc+1] + int2v[i+1][nc+1] + int2v[i][nc] + int2v[i-1][nc] + int2v[i+1][nc]) ; } /* End i loop */ for ( j=1 ; j<=nc ; j++ ) { int3v[0][j] = vstay*int2v[0][j] + ((1-vstay)/8)*(int2v[0][j-1] + int2v[0][j+1] + int2v[1][j-1] + int2v[1][j] + int2v[1][j+1]); int3v[nr+1][j] = vstay*int2v[nr+1][j] + ((1-vstay)/8)*(int2v[nr+1][j-1] + int2v[nr+1][j+1] + int2v[nr][j] + int2v[nr][j-1] + int2v[nr ][j+1]); /* That finishes the border interiors. Next, the plantation proper */ for ( i=1 ; i<=nr ; i++ ) { int3v[i][j] = vstay*int2v[i][j] + ((1-vstay)/8)*(int2v[i-1][j-1] + int2v[i-1][j] + int2v[i-1][j+1] + int2v[i][j-1] + int2v[i][j+1] + int2v[i+1][j-1] + int2v[i+1][j] + int2v[i+1][j+1]); } /* End i loop */ } /* End j loop */ /* Finally, the four corners of the border. */ int3v[0][0] = vstay*int2v[0][0] + ((1-vstay)/8)*(int2v[0][1] + int2v[1][1] + int2v[1][0]) ; int3v[0][nc+1] = vstay*int2v[0][nc+1] + ((1-vstay)/8)*(int2v[0][nc] + int2v[1][nc+1] + int2v[1][nc]) ; int3v[nr+1][0] = vstay*int2v[nr+1][0] + ((1-vstay)/8)*(int2v[nr][0] + int2v[nr+1][1] + int2v[nr][1]) ; int3v[nr+1][nc+1] = vstay*int2v[nr+1][nc+1] + ((1-vstay)/8)*(int2v[nr][nc+1] + int2v[nr+1][nc] + int2v[nr][nc]); /* * End of diffusion. Next, CONVECTION. Add wind vector to * coordinates of each point (i,j) and deposit there the volatiles * that were at (i,j). * int1v is the volatile matrix after convection. */ /* Initialize int1v matrix. */ for ( i=0 ; i<=nr+1 ; i++ ) { for ( j=0 ; j<=nc+1 ; j++ ) int1v[i][j] = 0.0 ; } /* End i loop */ /* Use the wind vector to move volatiles. */ for ( i=0 ; i<=nr+1 ; i++ ) { for ( j=0 ; j<=nc+1 ; j++ ) { k = floor(j + wind[0]) ; m = floor(i + wind[1]) ; /* Reminder -- wind(0) is east-west; wind(1) is north-south. */ if ( (k >= 0) && (k <= nr+1) && (m >= 0) && (m <= nc+1) ) int1v[k][m] = int3v[i][j] ; } /* End j loop */ } /* End i loop */ /* * Finally, add the CREATION of volatiles. Creation of volatiles is * proportional to crown size and occurs only at points in the * plantation, not at border points. */ for ( i=1 ; i<=nr ; i++ ) { for ( j=1 ; j<=nc ; j++ ) { int1v[i][j] = int1v[i][j] + create[i][j] ; /* N.B. Create = 0 if there is no tree (it's dead) */ } /* End j loop */ } /* End i loop */ /* * Now for the part of the loop that iterates until convergence * is achieved. */ cnsum = 0.0 ; for ( i=0 ; i<=nr+1 ; i++ ) { for ( j=0 ; j<=nc+1 ; j++ ) cnsum = cnsum + fabs( int1v[i][j] - oldvol[i][j] ); } /* End i loop */ err = cnsum ; n=n+1 ; if (n>100) /* OK - this just breaks the while loop */ break; for (i=0 ; i<=nr+1 ; i++) { for (j=0 ; j<=nc+1 ; j++) oldvol[i][j] = int1v[i][j] ; } /* End i loop */ } /* end of while err > tol */ for (i=0 ; i<=nr+1 ; i++) { for (j=0 ; j<=nc+1 ; j++) newvol[i][j] = oldvol[i][j] ; } /* End i loop */ /* * THE RECURSION FOR NUMBER OF ADULT WEEVILS (for all ndays in the day loop, * where 60 <= nday <= 345) * * N.B. -- Only non-ovipositors move; 'ovipositors' must * cease egg laying before they can move. Brood adults * are not added to the adults until the end of the season. */ for ( i=1 ; i<=nr ; i++ ) { for ( j=1 ; j<=nc ; j++ ) { if (wee_tree[i][j] == 0) continue; /* Do this computation only if the tree is alive */ /* For convenience, shorten the notation */ xr0 = r0[i][j] ; xt0 = t0[i][j] ; xlength = length[i][j]; xdiam = diam[i][j]; xh = newh[i][j]; xxx = xdiam*xlength ; factor = xh/(PI*xxx); if ( factor < 0.01) factor = 0.01 ; else if (factor > 100.0) factor = 100.0 ; /* * Adult MORTALITY was imposed in the computation for ovipositors; * Take out the OVIPOSITORS, since they don't move. * They will be added back in at the end of the day. */ newad[i][j] = newad[i][j] - newovi[i][j]; /* * Next, REPULSION of non-ovipositors * Note: Overall categories (disjoint) are ovipositor, * non-ovipositor, emerge. Within non-ovipositor, the * sub-categories are die, stay, repel, diffuse. */ xxx = RAD1*(1.0 + RAD2*xr0*xr0 + RAD3*xt0*xt0)*(1.0 + RAD4*factor); if(xxx > 1.0) xxx = 1.0; repel[i][j] = chance_check(newad[i][j], xxx); /* * Among other things, this computation ensures that weevils who move * to trees that have very short leaders are quickly repelled */ newad[i][j] = newad[i][j] - repel[i][j]; } /* End for j */ } /* End for i */ /* * Calculate the random DIFFUSION of moving weevils -- only those not * repelled can diffuse */ for (i=1; i<=nr; i++) { for (j=1; j<=nc; j++) { xxx = 0.0; if(adjacent[i][j] > 0) xxx = 1.0/adjacent[i][j]; randif[i][j] = chance_check(newad[i][j], xxx*(1.0-wstay)); /* Here, randif = # of weevils that move from tree i,j to EACH * (live) adjacent tree * * Don't need to worry if the 'repelling' tree is dead; in that * case, newad = 0, so there is no diffusion * Also, newad = 0 on the border, so no diffusion from the border */ } /* end j loop */ } /* end i loop */ /* Move the diffusing weevils */ for (i=1; i<=nr; i++) { for (j=1; j<=nc; j++) { if (wee_tree[i][j] == 1) /* Tree is alive */ { newad[i][j] = newad[i][j] - randif[i][j]*adjacent[i][j] + randif[i-1][j-1] + randif[i-1][j] + randif[i-1][j+1] + randif[i][j-1] + randif[i][j+1] + randif[i+1][j-1] + randif[i+1][j] + randif[i+1][j+1] ; if(newad[i][j] < 0) newad[i][j] = 0; } /* End if tree is alive */ } /* End for j */ } /* End for i */ /* * Now DISTRIBUTE the weevils that have * been repelled from their present tree */ for ( i=1 ; i<=nr ; i++ ) { for ( j=1 ; j<=nc ; j++ ) { if (wee_tree[i][j] == 1) /* Tree is alive */ { sum = 0.0 ; for ( k=0 ; k<=2 ; k++ ) { for ( m=0 ; m<=2 ; m++ ) { xxx = newvol[i-1+k][j-1+m] ; if ( xxx > 0.0) attract[k][m] = a0[i-1+k][j-1+m]*sqrt(xxx); else attract[k][m] = 0.0 ; /* Note -- This computation works because a0 is zero on the * border; also, a0 = 0 on dead trees */ sum = sum + attract[k][m] ; } /* End for m */ } /* End for k */ sum = sum - attract[1][1] ; attract[1][1] = 0.0 ; /* * The weevils in question have decided to LEAVE; they are not * attracted to the tree they are on now */ for ( k=0 ; k<=2 ; k++ ) { for ( m=0 ; m<=2 ; m++ ) { if (sum != 0.0) {num1 = chance_check(repel[i][j], attract[k][m]/sum); newad[i-1+k][j-1+m] = newad[i-1+k][j-1+m] + num1; repel[i][j] = repel[i][j] - num1; sum = sum - attract[k][m]; } /* End if sum not equal to zero */ } /* End for m */ } /* End for k */ } /* End if tree is alive */ } /* End for j */ } /* End for i */ /* Finally, add back in the OVIPOSITORS */ for (i=1; i<=nr; i++) { for (j=1; j<=nc; j++) { if (wee_tree[i][j] == 1) newad[i][j] = newad[i][j] + newovi[i][j]; } /* end j loop */ } /* end i loop */ /* End Recursion for Adult Weevils */ if (noweevil == 0) /* Here, we have weevils in the stand and want to compute and print average numbers of weevils, per live tree, in various categories, every 5 days */ { #ifdef HOOP if ((nday % 5) == 0) { avnonovi = 0.0; avovi = 0.0; avegg = 0.0; avlarv = 0.0; avpupa = 0.0; avyoungad = 0.0; avnewegg = 0.0; avemerge = 0.0; avlivebrood = 0.0; for (i=1; i<=nr; i++) { for (j=1; j<=nc; j++) { if (wee_tree[i][j] == 1) { avovi = avovi + newovi[i][j]; avnonovi = avnonovi + newad[i][j]; avnewegg = avnewegg + totnewegg[i][j]; avemerge = avemerge + totemerge[i][j]; avlivebrood = avlivebrood + livebrood[i][j]; /* youngad = number of juveniles with at least 737 degree-days * numegg = number of eggs, less than 143 degree-days * numlarv = number of larvae, between 143 and 568 degree-days * numpupa = number of pupae, between 569 and 736 degree-days */ youngad[i][j] = 0.0 ; numegg[i][j] = 0.0 ; numlarv[i][j] = 0.0 ; numpupa[i][j] = 0.0 ; for (n = 0; n<=198; n++) { xnum = newjuv[i][j][n] ; xdays = newdeg[i][j][n] ; if (xdays <= 142.0) numegg[i][j] = numegg[i][j] + xnum ; else if (xdays <= 568.0) numlarv[i][j] = numlarv[i][j] + xnum ; else if (xdays <= 736.0) numpupa[i][j] = numpupa[i][j] + xnum ; else youngad[i][j] = youngad[i][j] + xnum; } /* end for n */ avegg = avegg + numegg[i][j]; avlarv = avlarv + numlarv[i][j]; avpupa = avpupa + numpupa[i][j]; avyoungad = avyoungad + youngad[i][j]; } /* End if tree is alive */ } /* end j loop */ } /* end i loop */ /* Correct for the fact that newad includes both ovi and nonovi */ avnonovi = avnonovi - avovi ; if (numtrees > 0) { avovi = avovi/numtrees; avnonovi = avnonovi/numtrees; avnewegg = avnewegg/numtrees ; avemerge = avemerge/numtrees ; avlivebrood = avlivebrood/numtrees ; avegg = avegg/numtrees; avlarv = avlarv/numtrees; avpupa = avpupa/numtrees; avyoungad = avyoungad/numtrees; } else { avovi = 0.0; avnonovi = 0.0; avnewegg = 0.0; avemerge = 0.0; avlivebrood = 0.0; avegg = 0.0; avlarv = 0.0; avpupa = 0.0; avyoungad = 0.0; } /* Print summary statistics for the stand as a whole every 5 days. * Tree values: sumd, sumattk, percent, summerch, avhgt are zero before * nday = 365. During the season, all weevil averages are per live * tree, while after the season, variables avnewegg, avemerge and * avlivebrood are recomputed for trees whose leaders were 'attacked', * i.e., killed by weevils */ fprintf (errfile, "%4d %2d %3d %4d %2d", nsim, noweevil, numcln, pattern, weedist); fprintf (errfile, "%5.1f %5.1f %5.1f %2d", approxha, GRID, tmean, rancln); fprintf (errfile, "%3d %2d %5.2f %5.2f", repnum, izero, g0mean, R0MEAN); fprintf (errfile, "%5.2f %5.2f %3d %4d", T0MEAN, A0MEAN, nyear, nday); fprintf (errfile, "%7.1f %5d %5d %5.2f", totdeg, nr*nc, numtrees, sumd); fprintf (errfile, "%6d %6.1f %8.1f", sumattk, percent, summerch); fprintf (errfile, "%6.2f %6.2f %6.2f %6.2f", avhgt, avegg, avlarv, avpupa); fprintf (errfile, "%6.2f %6.2f %6.2f", avyoungad, avovi, avnonovi); fprintf (errfile, "%5.2f %5.1f %7.2f", A0VAR, multre1, avnewegg); fprintf (errfile, "%6.2f %6.2f \n", avemerge, avlivebrood); if (pattern != 0) { for (m=1; m<=numcln; m++) { avnonovi = 0.0; avovi = 0.0; avnewegg = 0.0; avemerge = 0.0; avlivebrood = 0.0; avegg = 0.0; avlarv = 0.0; avpupa = 0.0; avyoungad = 0.0; for (i=1; i<=nr; i++) { for (j=1; j<=nc; j++) { if (wee_tree[i][j] == 1 && treetype[i][j] == m) { avovi = avovi + newovi[i][j]; avnonovi = avnonovi + newad[i][j]; avnewegg = avnewegg + totnewegg[i][j]; avemerge = avemerge + totemerge[i][j]; avlivebrood = avlivebrood + livebrood[i][j]; youngad[i][j] = 0.0; numegg[i][j] = 0.0; numlarv[i][j] = 0.0; numpupa[i][j] = 0.0; for (n = 0; n<=198; n++) { xnum = newjuv[i][j][n]; xdays = newdeg[i][j][n]; if (xdays <= 142.0) numegg[i][j] = numegg[i][j] + xnum; else if (xdays <= 568.0) numlarv[i][j] = numlarv[i][j] + xnum; else if (xdays <= 736.0) numpupa[i][j] = numpupa[i][j] + xnum; else youngad[i][j] = youngad[i][j] + xnum; } /* end for n */ avegg = avegg + numegg[i][j]; avlarv = avlarv + numlarv[i][j]; avpupa = avpupa + numpupa[i][j]; avyoungad = avyoungad + youngad[i][j]; } /* End if tree is alive and clone m */ } /* End loop on j */ } /* End loop on i */ avnonovi = avnonovi - avovi; if (clntree[m] != 0) { avovi = avovi/clntree[m]; avnonovi = avnonovi/clntree[m]; avnewegg = avnewegg/clntree[m] ; avemerge = avemerge/clntree[m] ; avlivebrood = avlivebrood/clntree[m] ; avegg = avegg/clntree[m]; avlarv = avlarv/clntree[m]; avpupa = avpupa/clntree[m]; avyoungad = avyoungad/clntree[m]; } else { avovi = 0.0; avnonovi = 0.0; avnewegg = 0.0; avemerge = 0.0; avlivebrood = 0.0; avegg = 0.0; avlarv = 0.0; avpupa = 0.0; avyoungad = 0.0; } /* End if clntree[m] != 0 */ /* Print output re each clone every 5 days */ fprintf (errfile, "%4d %2d %3d %4d %2d", nsim, noweevil, numcln, pattern, weedist); fprintf (errfile, "%5.1f %5.1f %5.1f %2d", approxha, GRID, tmean, rancln); fprintf (errfile, "%3d %2d %5.2f %5.2f", repnum, m, g0cl[m], r0cl[m]); fprintf (errfile, "%5.2f %5.2f %3d %4d", t0cl[m], a0cl[m], nyear, nday); fprintf (errfile, "%7.1f %5d %5d", totdeg, startclntree[m], clntree[m]); fprintf (errfile, "%5.2f %6d %6.1f",sumdcl[m],sumattkcl[m],percentcl[m]); fprintf (errfile, "%8.1f %6.2f", summerchcl[m], avhgtcl[m]); fprintf (errfile, "%6.2f %6.2f %6.2f", avegg, avlarv, avpupa); fprintf (errfile, "%6.2f %6.2f %6.2f", avyoungad, avovi, avnonovi); fprintf (errfile, "%5.2f %5.1f", A0VAR, multre1); fprintf (errfile, "%7.2f %6.2f", avnewegg, avemerge); fprintf (errfile, "%6.2f \n", avlivebrood); } /* End loop on m */ } /* End if pattern != 0 */ } /* End if nday is divisible by 5 */ #endif } /* End if noweevil = 0 (i.e., there are weevils in the stand) */ } /* End of loop on nday (loop for the growing season) within each year (days 60-345) */ /* Assess the DAMAGE caused by weevils this year * 'Damage' = proportion of the leader diameter that has been * girdled by larvae in the feeding ring * * Then kill the leaders; the probability of leader death * is the 4th power of 'damage' (this formula fits the info * in the literature) * * In this code, weevil info is updated every day, so printed info is * always current, but printed tree info is current only after the * season ends (called nday = 365 here). All tree info printed * in the day loop above is for the preceding nday = 365 (from the * preceding year) */ nday = 365; for (i=1; i<=nr; i++) { for (j=1; j<=nc; j++) { if (wee_tree[i][j] == 1) /* tree is alive */ { leader[i][j] = 1 ; /* leader = 1 means the leader is alive; * now, kill those that are heavily damaged */ damage[i][j] = maxjuv[i][j]/(6.0*PI*diam[i][j]) ; if (damage[i][j] > 1.0) damage[i][j] = 1.0 ; xxx = damage[i][j] ; yyy = xxx*xxx*xxx*xxx ; urand = UniformRandom() ; /* urand is a uniform random number distributed U(0,1) */ if (urand <= yyy) leader[i][j] = 0 ; /* leader[i][j] = 0 if the leader of tree i,j is killed */ if (leader[i][j] == 0) { INT4 itn; itn = (i-1)*nc+(j-1)+1; if (tree_alive (itn) == TRUE) tree[itn]->data->itemp = 1; /* Tell TASS to mark tree's leader as killed */ } } /* End if tree is alive */ } /* end j loop */ } /* end i loop */ /* Compute the average damage per leader and the percentage of * surviving leaders in the stand as a whole and for each * clone, at the end of the season. * * Also compute total stand merchantable timber and average tree * height per live tree. * * Finally, compute this year's average total # of eggs laid, * average total # of emerged brood adults, and average * number of surviving brood adults, per ATTACKED tree, * ('attacked' = leader killed -- during the season, averages * were per LIVE tree) as of the end of the season (nday = 365) * These weevil values will all be zero if this is a 'no-weevil' * simulation. Also, if there are no weevils in the stand * this will be the only printout for the year (i.e., there * will be no printout during the season). */ percent = 0.0; sumd = 0.0; sumattk = 0; summerch = 0.0; avhgt = 0.0; avnewegg = 0.0; avemerge = 0.0; avlivebrood = 0.0; avovi = 0.0; avnonovi = 0.0; avegg = 0.0; avlarv = 0.0; avpupa = 0.0; avyoungad = 0.0; for (i=1; i<=nr; i++) { for (j=1; j<=nc; j++) { if (wee_tree[i][j] == 0) { wmerch[i][j] = 0.0; height[i][j] = 0.0; } else { itn = (i-1)*nc+(j-1)+1; merch (tree[itn]->data, &mrec); wmerch[i][j] = tree[itn]->data->vol_m; height[i][j] = tree[itn]->data->ht_t; percent = percent + leader[i][j] ; sumd = sumd + damage[i][j] ; summerch = summerch + wmerch[i][j] ; avhgt = avhgt + height[i][j] ; if (leader[i][j] == 0) { sumattk = sumattk + 1; avnewegg = avnewegg + totnewegg[i][j]; avemerge = avemerge + totemerge[i][j]; avlivebrood = avlivebrood + livebrood[i][j]; avnonovi = avnonovi + newad[i][j]; /* Assumes ovipositors switch to nonovi at end of season */ } } /* End else */ } /* end j loop */ } /* end i loop */ if (sumattk > 0) { /* Averages per attacked tree in the stand (this year) */ avnewegg = avnewegg/sumattk ; avemerge = avemerge/sumattk ; avlivebrood = avlivebrood/sumattk ; avnonovi = avnonovi/sumattk ; } if (numtrees > 0) { /* Averages per live tree in the stand (this year) */ percent = 100*percent/numtrees ; /* Percent of live trees with surviving leaders */ sumd = sumd/numtrees ; percentattk = 100.0*sumattk/numtrees ; /* Percent of live trees with killed leaders Ask KEN */ avhgt = avhgt/numtrees ; } else { percent = 0.0 ; sumd = 0.0 ; percentattk = 0.0 ; avhgt = 0.0 ; summerch = 0.0 ; } /* Print summary statistics for the stand as a whole. Note - some weevil averages are per attacked tree, and some are per live tree. See above. */ fprintf (errfile, "%4d %2d %3d %4d %2d", nsim, noweevil, numcln, pattern, weedist); fprintf (errfile, "%5.1f %5.1f %5.1f %2d", approxha, GRID, tmean, rancln); fprintf (errfile, "%3d %2d %5.2f %5.2f", repnum, izero, g0mean, R0MEAN); fprintf (errfile, "%5.2f %5.2f %3d %4d", T0MEAN, A0MEAN, nyear, nday); fprintf (errfile, "%7.1f %5d %5d %5.2f", totdeg, nr*nc, numtrees, sumd); fprintf (errfile, "%5.1f %5.1f %5.1f %5.1f", a0mult, g0mult, r0mult, t0mult); fprintf (errfile, "%6.1f %6.1f %8.1f", percentattk, percent, summerch); fprintf (errfile, "%6.2f %6.2f %6.2f %6.2f", avhgt, avegg, avlarv, avpupa); fprintf (errfile, "%6.2f %6.2f %6.2f", avyoungad, avovi, avnonovi); fprintf (errfile, "%5.2f %5.1f %7.2f", A0VAR, multre1, avnewegg); fprintf (errfile, "%6.2f %6.2f \n", avemerge, avlivebrood); if (pattern != 0) { for (m=1; m<=numcln; m++) { sumdcl[m] = 0.0 ; percentcl[m] = 0.0 ; summerchcl[m] = 0.0; avhgtcl[m] = 0.0; sumattkcl[m] = 0; /* # of trees of clone m with leaders killed this year */ avnewegg = 0.0; avemerge = 0.0; avlivebrood = 0.0; avovi = 0.0; avnonovi = 0.0; avegg = 0.0; avlarv = 0.0; avpupa = 0.0; avyoungad = 0.0; if (clntree[m] != 0) { for (i=1; i<=nr; i++) { for (j=1; j<=nc; j++) { if (wee_tree[i][j] == 1 && treetype[i][j] == m) { sumdcl[m] = sumdcl[m] + damage[i][j] ; percentcl[m] = percentcl[m] + leader[i][j] ; summerchcl[m] = summerchcl[m] + wmerch[i][j] ; avhgtcl[m] = avhgtcl[m] + height[i][j] ; if (leader[i][j] == 0) { sumattkcl[m] = sumattkcl[m] + 1; avnewegg = avnewegg + totnewegg[i][j]; avemerge = avemerge + totemerge[i][j]; avlivebrood = avlivebrood + livebrood[i][j]; avnonovi = avnonovi + newad[i][j]; } /* wmerch[i][j] and height[i][j] are zero if tree is dead */ } /* End if tree is alive and from clone m */ } /* End loop on j */ } /* End loop on i */ if (sumattkcl[m] > 0) { /* Averages per attacked tree of clone type m, in the stand (this year) */ avnewegg = avnewegg/sumattkcl[m] ; avemerge = avemerge/sumattkcl[m] ; avlivebrood = avlivebrood/sumattkcl[m] ; avnonovi = avnonovi/sumattkcl[m] ; } /* Averages per live tree of clone type m, in the stand */ sumdcl[m] = sumdcl[m]/clntree[m] ; percentattkcl[m] = 100*sumattkcl[m]/clntree[m] ; percentcl[m] = 100*percentcl[m]/clntree[m] ; avhgtcl[m] = avhgtcl[m]/clntree[m] ; } /* End if clntree[m] != 0 */ /* * Print summary statistics for each clone, if any. * Values of sumd, percentattk, percent, summerch, avhgt are zero * during the season */ fprintf (errfile, "%4d %2d %3d %4d %2d", nsim, noweevil, numcln, pattern, weedist); fprintf (errfile, "%5.1f %5.1f %5.1f %2d", approxha, GRID, tmean, rancln); fprintf (errfile, "%3d %2d %5.2f %5.2f", repnum, m, g0cl[m], r0cl[m]); fprintf (errfile, "%5.2f %5.2f %3d %4d", t0cl[m], a0cl[m], nyear, nday); fprintf (errfile, "%7.1f %5d %5d", totdeg, startclntree[m], clntree[m]); fprintf (errfile, "%5.1f %5.1f %5.1f %5.1f", a0mult, g0mult, r0mult, t0mult); fprintf (errfile, "%5.2f %6.1f %6.1f",sumdcl[m],percentattkcl[m],percentcl[m]); fprintf (errfile, "%8.1f %6.2f", summerchcl[m], avhgtcl[m]); fprintf (errfile, "%6.2f %6.2f %6.2f", avegg, avlarv, avpupa); fprintf (errfile, "%6.2f %6.2f %6.2f", avyoungad, avovi, avnonovi); fprintf (errfile, "%5.2f %5.1f", A0VAR, multre1); fprintf (errfile, "%7.2f %6.2f", avnewegg, avemerge); fprintf (errfile, "%6.2f \n", avlivebrood); } /* End loop on m */ } /* End if pattern not zero */ /* So far, this year's new brood have not participated in adult * activities. Now move new brood to newad. */ for (i=1; i<=nr; i++) { for (j=1; j<=nc; j++) { newad[i][j] = newad[i][j] + livebrood[i][j]; livebrood[i][j] = 0; } } /* Impose over-winter weevil mortality */ /* Compute the average number of degree-days after emergence */ totdeg = totdeg - 888.0; if (totdeg < 0.0) totdeg = 0.0 ; /* Functional form of wintersurv from Sullivan (1960). Value of DEG1 is chosen so, on average, wintersurv = 0.8 when average number of degree-days is 1700; i.e., when av temp = 10.75, as it is near Vancouver */ wintersurv = 1.0 - exp(-DEG1*totdeg); for (k=1; k<=nr; k++) { for (m=1; m<=nc; m++) newad[k][m] = chance_check(newad[k][m], wintersurv); } /* End k loop. */ /* Other weevil categories are set to zero at the beginning of each year */ #ifdef DEBUG printf ("weevils.c: weevil_year(): end\n"); #endif /* end of weevil_year() */ } /* End of loop on nyear */ void weevil_end () { /* * This function frees memory allocated at the start of the weevil * movement process. This function should be called once at the end * of the TASS run. */ short int i, j; for (i = 0; i < nr+2; i++) free (crown[i]); free (crown); for (i = 0; i < nr+2; i++) free (length[i]); free (length); for (i = 0; i < nr+2; i++) free (diam[i]); free (diam); for (i = 0; i < nr+2; i++) free (r0[i]); free (r0); for (i = 0; i < nr+2; i++) free (a0[i]); free (a0); for (i = 0; i < nr+2; i++) free (t0[i]); free (t0); for (i = 0; i < nr+2; i++) free (g0[i]); free (g0); for (i = 0; i < nr+2; i++) free (int1v[i]); free (int1v); for (i = 0; i < nr+2; i++) free (int2v[i]); free (int2v); for (i = 0; i < nr+2; i++) free (int3v[i]); free (int3v); for (i = 0; i < nr+2; i++) free (leader[i]); free (leader); for (i = 0; i < nr+2; i++) free (livebrood[i]); free (livebrood); for (i = 0; i < nr+2; i++) free (newad[i]); free (newad); for (i = 0; i < nr+2; i++) free (oldad[i]); free (oldad); for (i = 0; i < nr+2; i++) free (newovi[i]); free (newovi); for (i = 0; i < nr+2; i++) free (oldovi[i]); free (oldovi); for (i = 0; i < nr+2; i++) { for (j = 0; j < nc+2; j++) free (newjuv[i][j]); free (newjuv[i]); } free (newjuv); for (i = 0; i < nr+2; i++) { for (j = 0; j < nc+2; j++) free (oldjuv[i][j]); free (oldjuv[i]); } free (oldjuv); for (i = 0; i < nr+2; i++) free (maxjuv[i]); free (maxjuv); for (i = 0; i < nr+2; i++) free (emerge[i]); free (emerge); for (i = 0; i < nr+2; i++) { for (j = 0; j < nc+2; j++) free (newdeg[i][j]); free (newdeg[i]); } free (newdeg); for (i = 0; i < nr+2; i++) { for (j = 0; j < nc+2; j++) free (olddeg[i][j]); free (olddeg[i]); } free (olddeg); for (i = 0; i < nr+2; i++) free (newh[i]); free (newh); for (i = 0; i < nr+2; i++) free (oldh[i]); free (oldh); for (i = 0; i < nr+2; i++) free (neweh[i]); free (neweh); for (i = 0; i < nr+2; i++) free (oldeh[i]); free (oldeh); for (i = 0; i < nr+2; i++) free (newvol[i]); free (newvol); for (i = 0; i < nr+2; i++) free (oldvol[i]); free (oldvol); for (i = 0; i < nr+2; i++) free (create[i]); free (create); for (i = 0; i < nr+2; i++) free (repel[i]); free (repel); for (i = 0; i < nr+2; i++) free (postlarv2[i]); free (postlarv2); for (i = 0; i < nr+2; i++) free (youngad[i]); free (youngad); for (i = 0; i < nr+2; i++) free (numegg[i]); free (numegg); for (i = 0; i < nr+2; i++) free (numlarv[i]); free (numlarv); for (i = 0; i < nr+2; i++) free (numpupa[i]); free (numpupa); for (i = 0; i < nr+2; i++) free (totjuv[i]); free (totjuv); for (i = 0; i < nr+2; i++) free (randif[i]); free (randif); for (i = 0; i < nr+2; i++) free (damage[i]); free (damage); for (i = 0; i < nr+2; i++) free (wee_tree[i]); free (wee_tree); for (i = 0; i < nr+2; i++) free (adjacent[i]); free (adjacent); for (i = 0; i < nr+2; i++) free (treetype[i]); free (treetype); free (clntree); free (startclntree); free (startw); } /* End weevil_end */ // example usage: // num = chance_check (xnum, 1.0 - death); static short chance_check (short num_things, float chance_success) { short i; short num_success; num_success = 0; for (i = 0; i < num_things; i++) { if (UniformRandom () < chance_success) num_success++; } return num_success; } // example usage: // num = PoissonRandom (mu); /* static short PoissonRandom (float mu) { short i, poinum; float xxx, yyy, sum1; xxx = UniformRandom (); yyy = exp(-mu); sum1 = yyy; for (i = 0; i < 1000; i++) { if (xxx <= sum1) { poinum = i; break; } else { yyy = yyy*mu/(i+1); sum1 = sum1 + yyy; } } return poinum; } */