/*** NAME CodPGR2 PURPOSE Module computes the cod population growth rate as a function of a parameter, assuming that resources are constant. An important assumption is that starvation of cod does not occur, i.e. the net-energy production never becomes negative Last modification: FHS - June 4, 2020 ***/ #include "globals.h" //======= PROGRAM SETTINGS ========================================================================================================= #define SOLVEFOR PGR //#define SOLVEFOR C_TSH #if(SOLVEFOR == C_TSH) #pragma message("Solving for: C_TSH") #elif(SOLVEFOR == PGR) #pragma message("Solving for: PGR") #endif /* *==================================================================================================================================== * DEFINITION OF PROBLEM DIMENSIONS AND NUMERICAL SETTINGS *==================================================================================================================================== */ // Dimension settings: Required #define EQUATIONS_DIM 1 #define EXTRAOUTPUT_DIM 7 #define PARAMETER_NR 35 // Numerical settings: Optional (default values adopted otherwise) #define ODESOLVE_MAX_STEP 1.0 // Largest step size in odesolver #define DYTOL 1.0E-7 // Variable tolerance #define RHSTOL 1.0E-8 // Function tolerance #define ALLOWNEGATIVE 1 // Negative solution components allowed? // Numerical settings #define JACOBIAN_STEP 1.0E-4 #define MIN_SURVIVAL 1.0E-9 // Survival at which individual is considered dead #define ZERO 1.0E-4 // Threshold for indentity with 0 /* *==================================================================================================================================== * DEFINITION OF ALIASES *==================================================================================================================================== */ // Define aliases for the parameters #define PGR parameter[ 0] // Population growth rate #define FL0 parameter[ 1] // 1st Feeding level coefficient #define FL1 parameter[ 2] // 2nd feeding level coef. #define FLL1 parameter[ 3] // Length at start switch FL #define FLL2 parameter[ 4] // Length at 50% switch FL #define YEAR parameter[ 5] #define C_SPAWNSET parameter[ 6] #define C_LWC parameter[ 7] // Length-weight scaling constant #define C_LWE parameter[ 8] // Length-weight scaling exponent #define C_REPROEFF parameter[ 9] // Gonad-offspring conversion efficiency #define C_BIRTHLENGTH parameter[10] // Length at birth #define C_MULENGTH parameter[11] // Size-dep. mortality characteristic size #define C_MATURELEN parameter[12] // Maturation length #define C_AGEFEEDING parameter[13] // Age at first feeding #define C_QJ parameter[14] // Juvenile condition target #define C_QA parameter[15] // Adult condition target #define C_DIGTIMEC parameter[16] // Digestion time scaling constant #define C_DIGTIMEE parameter[17] // Digestion time scaling exponent #define C_MAINTC parameter[18] // Maintenance scaling constant #define C_MAINTE parameter[19] // Maintenance scaling exponent #define C_CONVEFF parameter[20] // Conversion efficiency #define C_EGGMORT parameter[21] // Egg mortality #define C_MUSDC parameter[22] // Size-dependent mortality constant #define C_MUB parameter[23] // Background mortality #define C_FSTARTL parameter[24] // Size of start fishing vulnerability #define C_FHALFL parameter[25] // Size of 50% fishing vulnerability #define C_FDECL parameter[26] // Size of start decrease fishing vulnerability #define C_FDECHL parameter[27] // Size of 50% of decrease fishing vulnerability #define C_MINF parameter[28] // Minimum Fisheries retention large sizes #define C_FMORT parameter[29] // Fishing mortality #define C_Q10FACTOR parameter[30] // Q10 factor for cod energetics //additional parameters to manipulate the effect of sound disturbances #define C_MUS parameter[31] // Mortality multiplier #define C_TSH parameter[32] // Time spend hiding (proportion) #define C_EVI parameter[33] // Egg viability multiplier #define C_MSM parameter[34] // Maintenance multiplier /* *==================================================================================================================================== * DEFINITION OF NAMES AND DEFAULT VALUES OF THE PARAMETERS *==================================================================================================================================== */ // At least two parameters should be specified in this array char *parameternames[PARAMETER_NR] = {"PGR", "Flo0", "Flo1", "FlL1", "FlL2", "Year", "Spawn", "LambdaC", "LambdaE", "ThetaP", "L0", "Lsd", "Lm", "Aj", "qj", "qa", "XiC", "XiE", "RhoC", "RhoE", "ThetaA", "MuE", "MuJ", "Mu0", "Lh", "Lh50", "Ld", "Ld50", "Retm", "Muh", "Q10", "MuS", "Tsh", "Evi", "Msm"}; // These are the default parameters values double parameter[PARAMETER_NR] = {0.0, 1.0, 0.0, 0.39, 30.0, 250.0, 200.0, 0.01, 3.0, 0.5, 0.39, 3.68, 62.0, 22.0, 0.7, 1.2, 270.651 , -2.389, 0.022, 0.9124, 0.6, 0.03, 0.03, 0.003, 10.0, 34.0,58.0, 78.0, 0.55, 0.00124, 1.0, 0.0, 0.0, 0.0, 0.0}; /* *==================================================================================================================================== * DEFINITION OF THE SYSTEM OF EQUATIONS TO SOLVE *==================================================================================================================================== * * Specify the threshold determining the end point of each discrete life * stage in individual life history as function of i-state and environment * variables for all individuals of all populations in every life stage * * Notice that the first index of the variable 'istate[][]' refers to the * number of the structured population, the second index refers to the * number of the individual state variable. The interpretation of the latter * is up to the user. *==================================================================================================================================== */ double Sigmoid(double val, double low, double half) { double result = 0.0, nx; // Scale the dependent range between 0 and 3 nx = 1.5*(val - low)/(half - low); if (nx <= 0.0) result = 0.0; else if (nx <= 1.0) result = nx*nx*nx/6.0; else if (nx <= 2.0) result = -3.0*nx/2.0 + 3.0*nx*nx/2.0 - nx*nx*nx/3.0 + 0.5; else if (nx <= 3.0) result = 9.0*nx/2.0 - 3.0*nx*nx/2.0 + nx*nx*nx/6.0 - 3.5; else result = 1.0; return result; } /*==================================================================================================================================*/ // The ODE system defining the change in state variables during the growth period #define ODEDIM 7 static double StoredVals[ODEDIM]; static int OdeDim = ODEDIM; static double TotalMortLimit, R0, LifeRep, Fecundity; static int Mature, Spawning; #define TOTALMORT argument[0] #define AGE argument[1] #define LENGTH argument[2] #define RESERVES argument[3] #define GONADS argument[4] #define HARVESTINGYIELD argument[5] #define REALMORT argument[6] void WithinSeason(double t, double *argument, double *derivative) { double mass, bone, weight, fatratio; double htime, feedl; double maint, net_energy, q, kappa; double survival, mort, hmort; mass = C_LWC*pow(LENGTH, C_LWE); bone = mass/(1 + C_QJ); weight = bone + max(RESERVES + GONADS, 0.0); mort = PGR; mort += (1 + C_MUS)*C_MUB; // Size-independent background and sound mortality mort += (C_MUSDC > 0)*C_MUSDC*exp(-pow((LENGTH/C_MULENGTH),C_LWE)); // Size-dependent mortality hmort = (C_FMORT > 0)*C_FMORT*(Sigmoid(LENGTH, C_FSTARTL, C_FHALFL) - (1 - C_MINF)*Sigmoid(LENGTH, C_FDECL, C_FDECHL)); // Harvesting mortality survival = exp(-REALMORT); // use the survival without the PGR for calculations over the life cycle. PGR is not a real source of mortality. htime = C_DIGTIMEC*pow(LENGTH, C_DIGTIMEE); feedl = FL0 + FL1*Sigmoid(LENGTH,FLL1,FLL2); // The feeding level represents a simple T2FR; the value ranges between 0-1. feedl *= (1 - C_TSH); maint = C_Q10FACTOR*C_MAINTC*pow(weight, C_MAINTE); maint *= (1 + C_MSM); net_energy = C_CONVEFF*C_Q10FACTOR*(feedl / htime) - maint; if (Mature) q = C_QA; else q = C_QJ; fatratio = (RESERVES + GONADS)/bone; if (fatratio > q) kappa = 1.0 / (q + 1.0); else // Default recovery formula changed to a quadratic dependence on fatratio (below), such that recovery in // reversible mass is faster kappa = fatratio*fatratio / ((1.0 + q)*q*q); // dBones/da = dBones/dL * dL/da => dL/da = (dBones/da)/(dBones/dL) // Bones = (c*L^d)/(1 + qj) => dBones/dL = (c*d*L^(d-1))/(1 + qj) = d*Bones/L derivative[0] = mort + hmort; derivative[1] = 1.0; derivative[2] = kappa*net_energy*LENGTH/(C_LWE*bone); if (!Spawning) { derivative[3] = (1.0 - kappa)*net_energy; derivative[4] = 0.0; } else { derivative[3] = 0.0; derivative[4] = (1.0 - kappa)*net_energy; } derivative[5] = hmort*weight*survival; // use the survival without the PGR for the calculation of the yield. PGR is not a real source of mortality. derivative[6] = hmort + (mort - PGR); return; } double CheckAdultlen(double age, double *argument) { return max((LENGTH - C_MATURELEN), (REALMORT - TotalMortLimit)); } double CheckMinSurvival(double age, double *argument) { return (REALMORT - TotalMortLimit); } void PrintState(double *argument) { double vec[ODEDIM + 4]; //output in test run vec[0] = AGE; //age vec[1] = exp(-REALMORT); //survival memcpy(vec+2, argument + 2, (OdeDim - 2)*sizeof(double)); //all ODEs from length vec[ODEDIM] = R0; //R0 vec[ODEDIM + 1] = LifeRep; //LifeRep vec[ODEDIM + 2] = C_QJ*C_LWC*pow(LENGTH, C_LWE)/(1.0 + C_QJ); //condition vec[ODEDIM + 3] = Fecundity; //egg production PrettyPrintArray(NULL, ODEDIM + 4, vec); return; } /*==================================================================================================================================*/ // Routine specifying the system of equalities from which to solve for // the unknown variables int Equations(double *unknowns, double *result) { int years, odereturn; double tval, tnext, tend, argument[ODEDIM], bweight, timeinyear, bone; SOLVEFOR = unknowns[0]; //================================================================================================================================== // Set the initial point for the ODEs and flags memset(argument, 0, ODEDIM*sizeof(double)); TOTALMORT = (C_EGGMORT + PGR)*C_AGEFEEDING; AGE = C_AGEFEEDING; bweight = C_LWC*pow(C_BIRTHLENGTH, C_LWE); LENGTH = C_BIRTHLENGTH; RESERVES = (C_QJ) / (1.0 + C_QJ)*bweight; GONADS = 0.0; HARVESTINGYIELD = 0.0; REALMORT = C_EGGMORT*C_AGEFEEDING; TotalMortLimit = -log(MIN_SURVIVAL); Mature = 0; Spawning = 0; R0 = 0.0; LifeRep = 0.0; Fecundity = 0.0; tval = C_AGEFEEDING; tend = C_AGEFEEDING + (TotalMortLimit - (TOTALMORT - PGR)) / (C_MUB); // adjusted to calculate the full life history (without PGR related mortality term --- too slow??) tnext = C_SPAWNSET; if (TestRun) { STDOUT("\n\n Age Survival Length Reserves Gonads Yield Total_mort R0 LifeRep qj_bone Fecundity \n\n"); PrintState(argument); } while (tval < tend) { if (!Mature) // Juvenile fish { // Stopping at the end of a spawning period or year is not necessary, as the // seasonal cycle does not affect juveniles odereturn = odesolve(argument, OdeDim, &tval, tend, WithinSeason, CheckAdultlen); if (odereturn == FAILURE) { ErrorMsg(__FILE__, __LINE__, "Integration failed!"); return FAILURE; } if (TestRun) PrintState(argument); // Integration stops only through reaching maturation length or survival limit Mature = 1; // Determine in which year we are and whether a spawning setting moment is still coming up this year // If not set the year value 1 higher, because spawning is only going to happen in the next year years = (int)floor((tval + ZERO)/YEAR); if (tval >= (years*YEAR + C_SPAWNSET)) years++; tnext = (years*YEAR + C_SPAWNSET); // If minimum survival threshold is reaching stop the life history integration if (fabs(REALMORT - TotalMortLimit) < ZERO) break; } else // Mature fish { Fecundity = 0.0; odereturn = odesolve(argument, OdeDim, &tval, tnext, WithinSeason, CheckMinSurvival); if (odereturn == FAILURE) { ErrorMsg(__FILE__, __LINE__, "Integration failed!"); return FAILURE; } if (TestRun) PrintState(argument); if (odereturn == SUCCES_ODE_MAXTIME) // Reached time of spawn set or end of year { timeinyear = tval - years*YEAR; if (Spawning) // Must be end of year { if (fabs(tval - (years + 1)*YEAR) > ZERO) { ErrorMsg(__FILE__, __LINE__, "Integration failed to stop at the end of spawning period: T = %12.6f (should be %12.6f!)", tval, (years + 1)*YEAR); PrintState(argument); return FAILURE; } Fecundity = max((1.0 - C_EVI)*(C_REPROEFF*GONADS / bweight), 0.0); // Add the offspring produced to R0 R0 += Fecundity*exp(-TOTALMORT); //LifeRep += exp(-REALMORT); LifeRep += Fecundity*exp(-REALMORT); //*exp(-REALMORT); GONADS = 0.0; Spawning = 0; years++; tnext = years*YEAR + C_SPAWNSET; if (TestRun) PrintState(argument); } else { if ((fabs(timeinyear - C_SPAWNSET) > ZERO)) { ErrorMsg(__FILE__, __LINE__, "Integration failed to stop at the start of spawning period: T = %12.6f (should be %12.6f!, tnext = %12.6f)", tval, years*YEAR + C_SPAWNSET, tnext); PrintState(argument); return FAILURE; } bone = C_LWC*pow(LENGTH, C_LWE) / (1.0 + C_QJ); //spawn setting depends on the body condition if(RESERVES < C_QJ*bone) { Spawning = 0; //spawning will not occur this year, so we need to step to the setspawning event next year years++; tnext = years*YEAR + C_SPAWNSET; } else { Spawning = 1; //part of the mass in the reserves is transferred to the gonads GONADS = RESERVES - C_QJ*bone; RESERVES = C_QJ*bone; //spawning will occur this year, so we step to the spawning event this year tnext = (years + 1)*YEAR; } if (TestRun) PrintState(argument); } } // If minimum survival threshold is reaching stop the life history integration if (fabs(REALMORT - TotalMortLimit) < ZERO) break; } } if (TestRun) PrintState(argument); memcpy(StoredVals, argument, ODEDIM*sizeof(double)); //================================================================================================================================== // Compute the final values of the fixed point equation F(y)=0, if (result) result[0] = (R0 - 1.0); return SUCCES; } /*==================================================================================================================================*/ // Define all variables to be written to the output file (column-organized ASCII file) int DefineExtraOutput(double *argument, double *ExtraOutput) { // Invoke the routine that sets the right-hand side for setting the output variables. if (Equations(argument, NULL) == FAILURE) return FAILURE; ExtraOutput[0] = StoredVals[ 1]; // Age ExtraOutput[1] = exp(-StoredVals[0]); // Survival ExtraOutput[2] = StoredVals[ 2]; // Length ExtraOutput[3] = StoredVals[ 3]; // Reserves ExtraOutput[4] = StoredVals[ 4]; // Gonads ExtraOutput[5] = StoredVals[ 5]; // Harvesting yield ExtraOutput[6] = R0; ExtraOutput[7] = LifeRep; // LifeTimeReproductiveOutput (based on real mort) return SUCCES; } /* *==================================================================================================================================== * Default values for several indices used in the bifurcation: * * Bifparone : Index of the first free parameter in the bifurcation * Bifpartwo : Index of the second free parameter in the bifurcation * EXTfun : Index of the element in the ExtraOutput[] vector, specifying * the function value for which to monitor or continue a maximum * or minimum value along the curve * EXTpar : Index of the parameter, with respect to which to determine the * derivative of the function value contained in ExtraOutput[EXTfun] * along the curve * * These indices can also be set via command-line options *==================================================================================================================================== */ void SetBifPars(void) { Bifparone = 1; Bifpartwo = 2; return; } /*==================================================================================================================================*/