Commit 326a3b75 authored by Olivier Bonnefon's avatar Olivier Bonnefon
Browse files

AEDES interAnual model update

parent 01b2aaf7
...@@ -40,6 +40,7 @@ alphaw1=1; ...@@ -40,6 +40,7 @@ alphaw1=1;
//h0/364~0.1 //h0/364~0.1
//h0=36.4 //h0=36.4
//real lifeexpectancy=1.0/12.0; //real lifeexpectancy=1.0/12.0;
real h0=36.4; real h0=36.4;
//deathrate=1/lifeexpected //deathrate=1/lifeexpected
real deathrate=12; real deathrate=12;
...@@ -61,7 +62,9 @@ if (scenario==scenarioHight){ ...@@ -61,7 +62,9 @@ if (scenario==scenarioHight){
} }
real K=1.0; real K=1.0;
real Ru2=(Rstar-deathrate)/K; //real Ru2=(Rstar-deathrate)/K;
real maxGt=0.11*365.0;
real Ru2=maxGt - deathrate;
include "defOutputGen.edp"; include "defOutputGen.edp";
real[int] u6(6); real[int] u6(6);
//real[int] grad6(6); //real[int] grad6(6);
......
...@@ -99,39 +99,39 @@ if (numDay > 59){ ...@@ -99,39 +99,39 @@ if (numDay > 59){
LoadVec(GT14[],DATASAFRAN +"GT_BY_DAY/"+string(CURYEAR)+"/mesh14s"+string(numDay)); LoadVec(GT14[],DATASAFRAN +"GT_BY_DAY/"+string(CURYEAR)+"/mesh14s"+string(numDay));
//update R //update R
RR0=Rstar*GT0*GP0; RR0=GT0*GP0;
RR1=Rstar*GT1*GP1; RR1=GT1*GP1;
RR2=Rstar*GT2*GP2; RR2=GT2*GP2;
RR3=Rstar*GT3*GP3; RR3=GT3*GP3;
RR4=Rstar*GT4*GP4; RR4=GT4*GP4;
RR5=Rstar*GT5*GP5; RR5=GT5*GP5;
RR6=Rstar*GT6*GP6; RR6=GT6*GP6;
RR7=Rstar*GT7*GP7; RR7=GT7*GP7;
RR8=Rstar*GT8*GP8; RR8=GT8*GP8;
RR9=Rstar*GT9*GP9; RR9=GT9*GP9;
RR10=Rstar*GT10*GP10; RR10=GT10*GP10;
RR11=Rstar*GT11*GP11; RR11=GT11*GP11;
RR12=Rstar*GT12*GP12; RR12=GT12*GP12;
RR13=Rstar*GT13*GP13; RR13=GT13*GP13;
RR14=Rstar*GT14*GP14; RR14=GT14*GP14;
//if (numDay%10==0) //if (numDay%10==0)
// plot(RR0,RR1,RR2,RR3,RR4,RR5,RR6,RR7,RR8,RR9,RR10,RR11,RR12,RR13,RR14,cmm="rr "+string(numDay),value=true,fill=true,wait=true); // plot(RR0,RR1,RR2,RR3,RR4,RR5,RR6,RR7,RR8,RR9,RR10,RR11,RR12,RR13,RR14,cmm="rr "+string(numDay),value=true,fill=true,wait=true);
//update gammaw1s1 //update gammaw1s1
gammaw1s0=Rstar*GT0*(1-GP0); gammaw1s0=GT0*(1-GP0);
gammaw1s1=Rstar*GT1*(1-GP1); gammaw1s1=GT1*(1-GP1);
gammaw1s2=Rstar*GT2*(1-GP2); gammaw1s2=GT2*(1-GP2);
gammaw1s3=Rstar*GT3*(1-GP3); gammaw1s3=GT3*(1-GP3);
gammaw1s4=Rstar*GT4*(1-GP4); gammaw1s4=GT4*(1-GP4);
gammaw1s5=Rstar*GT5*(1-GP5); gammaw1s5=GT5*(1-GP5);
gammaw1s6=Rstar*GT6*(1-GP6); gammaw1s6=GT6*(1-GP6);
gammaw1s7=Rstar*GT7*(1-GP7); gammaw1s7=GT7*(1-GP7);
gammaw1s8=Rstar*GT8*(1-GP8); gammaw1s8=GT8*(1-GP8);
gammaw1s9=Rstar*GT9*(1-GP9); gammaw1s9=GT9*(1-GP9);
gammaw1s10=Rstar*GT10*(1-GP10); gammaw1s10=GT10*(1-GP10);
gammaw1s11=Rstar*GT11*(1-GP11); gammaw1s11=GT11*(1-GP11);
gammaw1s12=Rstar*GT12*(1-GP12); gammaw1s12=GT12*(1-GP12);
gammaw1s13=Rstar*GT13*(1-GP13); gammaw1s13=GT13*(1-GP13);
gammaw1s14=Rstar*GT14*(1-GP14); gammaw1s14=GT14*(1-GP14);
CURYEAR=CURYEARSAV; CURYEAR=CURYEARSAV;
...@@ -4,7 +4,7 @@ ...@@ -4,7 +4,7 @@
#include <fstream> #include <fstream>
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include <math.h>
using namespace std; using namespace std;
#include "../defRasterSafran.h" #include "../defRasterSafran.h"
...@@ -230,7 +230,9 @@ void computeT0usingGDD(double T0,double GDD_Threshold,double *pDataT,int dim, do ...@@ -230,7 +230,9 @@ void computeT0usingGDD(double T0,double GDD_Threshold,double *pDataT,int dim, do
} }
} }
/*
*
*previous version piecewise linear
void computeGT(double *pDataT,int dim,double *pGT){ void computeGT(double *pDataT,int dim,double *pGT){
double *pCurGT=pGT; double *pCurGT=pGT;
double *pAux; double *pAux;
...@@ -255,8 +257,53 @@ void computeGT(double *pDataT,int dim,double *pGT){ ...@@ -255,8 +257,53 @@ void computeGT(double *pDataT,int dim,double *pGT){
(*(pCurGT++))/=24.0; (*(pCurGT++))/=24.0;
} }
} }
*/
/*
10.804361785896965 0.080201779776833 12.438168144101718
Columns 4 through 5
40.439255430520049 -0.145790284196338
*/
double psi=10.804361785896965;
double alpha=0.080201779776833;
double delta=12.438168144101718;
double Tmax=40.439255430520049;
double lambda= -0.145790284196338 ;
double rhofct(double T){
double res=psi*(exp(alpha*T)-exp(alpha*Tmax-(Tmax-T)/delta))+lambda;
if (res<0)
res=0;
return res;
}
void computeGT(double *pDataT,int dim,double *pGT){
double *pCurGT=pGT;
double *pAux;
for(int j=0;j<365;j++){
pAux=pCurGT;
for (long int i=0;i<dim;i++)
*(pAux++)=0;
for (int h=0;h<24;h++){
pAux=pCurGT;
for (long int i=0;i<dim;i++){
double temp=*(pDataT++);
*(pAux++)+=temp;
}
}
for (long int i=0;i<dim;i++){
(*pCurGT)/=24.0;
*pCurGT=rhofct(*pCurGT);
pCurGT++;
}
}
}
void computeAverageT(double Period,double *pDataT,int dim,double *pDataAverage){ void computeAverageT(double Period,double *pDataT,int dim,double *pDataAverage){
double *pBuf=(double *)calloc(dim,sizeof(double)); double *pBuf=(double *)calloc(dim,sizeof(double));
double *pAux=pBuf; double *pAux=pBuf;
...@@ -367,7 +414,7 @@ int main(){ ...@@ -367,7 +414,7 @@ int main(){
readTemperatures(pDatas,meshX,meshY,dim,numA); readTemperatures(pDatas,meshX,meshY,dim,numA);
computeGDD(10,pDatas,dim, pT0J); /*computeGDD(10,pDatas,dim, pT0J);
saveField(pT0J,numA,"GDD_BY_DAY",365-59); saveField(pT0J,numA,"GDD_BY_DAY",365-59);
seuilGDD(294,dim,pT0J,pT0JSeuil); seuilGDD(294,dim,pT0J,pT0JSeuil);
saveField(pT0JSeuil,numA,"GDD_294",365-59); saveField(pT0JSeuil,numA,"GDD_294",365-59);
...@@ -376,7 +423,7 @@ int main(){ ...@@ -376,7 +423,7 @@ int main(){
seuilGDD(246,dim,pT0J,pT0JSeuil); seuilGDD(246,dim,pT0J,pT0JSeuil);
saveField(pT0JSeuil,numA,"GDD_246",365-59); saveField(pT0JSeuil,numA,"GDD_246",365-59);
continue; continue;*/
computeGT(pDatas,dim,pGT); computeGT(pDatas,dim,pGT);
saveField(pGT,numA,"GT_BY_DAY",365); saveField(pGT,numA,"GT_BY_DAY",365);
// computeT0usingGDD(10,5,pDatas,dim, pT0J); // computeT0usingGDD(10,5,pDatas,dim, pT0J);
......
...@@ -26,12 +26,8 @@ string FIELD; ...@@ -26,12 +26,8 @@ string FIELD;
//FIELD = "END_DIAPOSE_BY_DAY/"; //FIELD = "END_DIAPOSE_BY_DAY/";
//FIELD = "BEGIN_DIAPOSE_BY_DAY/"; //FIELD = "BEGIN_DIAPOSE_BY_DAY/";
//FIELD = "GDDTHRESHOLD_BY_DAY/"; //FIELD = "GDDTHRESHOLD_BY_DAY/";
FIELD = "GDD_BY_DAY/"; FIELD = "GT_BY_DAY/";
FIELD = "GDD_294/"; //FIELD = "GP_BY_DAY/";
FIELD = "GDD_269/";
FIELD = "GDD_246/";
//FIELD = "GT_BY_DAY/";
FIELD = "GP_BY_DAY/";
int step=1; int step=1;
// LoadVec(T[],"T_BY_HEURE/meshT0"); // LoadVec(T[],"T_BY_HEURE/meshT0");
// plot(T,wait=true,cmm="T0"); // plot(T,wait=true,cmm="T0");
...@@ -68,7 +64,7 @@ if (FIELD == "PHOTO_PERIOD_BY_DAY/"){ ...@@ -68,7 +64,7 @@ if (FIELD == "PHOTO_PERIOD_BY_DAY/"){
if (FIELD == "GP_BY_DAY/"){ if (FIELD == "GP_BY_DAY/"){
for (int j=210;j<365;j++){ for (int j=1;j<365;j++){
LoadVec(T0[],DIRNAME +FIELD+"/gp0s"+string(j)); LoadVec(T0[],DIRNAME +FIELD+"/gp0s"+string(j));
LoadVec(T1[],DIRNAME +FIELD+"/gp1s"+string(j)); LoadVec(T1[],DIRNAME +FIELD+"/gp1s"+string(j));
LoadVec(T2[],DIRNAME +FIELD+"/gp2s"+string(j)); LoadVec(T2[],DIRNAME +FIELD+"/gp2s"+string(j));
...@@ -111,8 +107,7 @@ if (FIELD == "GT_BY_DAY/"){ ...@@ -111,8 +107,7 @@ if (FIELD == "GT_BY_DAY/"){
} }
if (FIELD == "GDDTHRESHOLD_BY_DAY/" || FIELD == "GDD_BY_DAY/" || if (FIELD == "GDDTHRESHOLD_BY_DAY/"){
FIELD== "GDD_294/" ||FIELD== "GDD_269/" ||FIELD== "GDD_246/" ){
for (int j=1;j<365-59;j++){ for (int j=1;j<365-59;j++){
LoadVec(T0[],DIRNAME +FIELD+string(numA)+"/mesh0s"+string(j)); LoadVec(T0[],DIRNAME +FIELD+string(numA)+"/mesh0s"+string(j));
LoadVec(T1[],DIRNAME +FIELD+string(numA)+"/mesh1s"+string(j)); LoadVec(T1[],DIRNAME +FIELD+string(numA)+"/mesh1s"+string(j));
...@@ -129,7 +124,7 @@ if (FIELD == "GDDTHRESHOLD_BY_DAY/" || FIELD == "GDD_BY_DAY/" || ...@@ -129,7 +124,7 @@ if (FIELD == "GDDTHRESHOLD_BY_DAY/" || FIELD == "GDD_BY_DAY/" ||
LoadVec(T12[],DIRNAME +FIELD+string(numA)+"/mesh12s"+string(j)); LoadVec(T12[],DIRNAME +FIELD+string(numA)+"/mesh12s"+string(j));
LoadVec(T13[],DIRNAME +FIELD+string(numA)+"/mesh13s"+string(j)); LoadVec(T13[],DIRNAME +FIELD+string(numA)+"/mesh13s"+string(j));
LoadVec(T14[],DIRNAME +FIELD+string(numA)+"/mesh14s"+string(j)); LoadVec(T14[],DIRNAME +FIELD+string(numA)+"/mesh14s"+string(j));
plot(T0,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,cmm=FIELD+" "+string(j),value=true,fill=true,wait=false); plot(T0,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,cmm="GDDTHRESHOLD_BY_DAY j "+string(j),value=true,fill=true,wait=false);
} }
} }
......
...@@ -26,6 +26,8 @@ Ce repertoire contient les scripts necessaires à l'exploitation des données sa ...@@ -26,6 +26,8 @@ Ce repertoire contient les scripts necessaires à l'exploitation des données sa
3.2) Cas d'un seul maillage representant la france (ie sous repertoire FRANCE) 3.2) Cas d'un seul maillage representant la france (ie sous repertoire FRANCE)
//fichier computeIjxy.c
permet de tracer le nuage de point (Ijxy,Tj)
......
#include "stdio.h"
#include "string.h"
#include <stdlib.h>
#include <sys/dir.h>
#include <sys/stat.h>
#include <sys/types.h>
#include "math.h"
#define NB_PT_SAFRAN 9892
#define NB_TIME_SAFRAN 8760
#define FIRST_YEARS 2011
#define LAST_YEARS 2012
#define T_ZERO_OBSOLU -273.15
#define PATH_TO_DATA_SAFRAN "/home/olivierb/solvers/trunk/SandBox/FF/AedesAlbopictus/intraAnnuel/donnees_SAFRAN/AUTOROUTES/DATA/"
//11.128007398517115 0.078865315431034 12.651026892598505 42.345187120388388 -0.145204446756075
double psi=11.128007398517115;
double alpha=0.078865315431034;
double delta=12.651026892598505;
double Tmax=42.345187120388388;
double lambda=-0.145204446756075 ;
//interval de temps:
/*
on veut les dates comprise entre le 1 mars et le 31 dec
Les donnees de 2011 commence le 1aout 2011.
on garde les dates du 1aout au 31 dec : de 0 a 3666
on ignore les dates du 1j au 1 mars ie de 3666 a 59*24+3666=5082
*/
double Tmin=0;
double Vmin;
double Vmax=-100;
double rhofct(double T){
double res=psi*(exp(alpha*T)-exp(alpha*Tmax-(Tmax-T)/delta))+lambda;
if (T<Tmin){
Tmin=T;
Vmin=res;
}
if (Vmax<res)
Vmax=res;
return res;
}
#define PERIOD_HEURE 24
void safranToIxyTday( ){
char inputDataSafran[512];
double *pValT=(double*)malloc(NB_PT_SAFRAN*sizeof(double));
int step=1;
double * pTday=(double *)calloc(NB_PT_SAFRAN,sizeof(double));
double * pIjxy=(double *)calloc(NB_PT_SAFRAN,sizeof(double));
FILE * foutput=fopen("IjxyTj.txt","w");
int cmp=0;
for (int numA=FIRST_YEARS;numA<LAST_YEARS;numA++){
int k;
for (int step=1; step<NB_TIME_SAFRAN+1; step++){
if (step > 3666 && step < 5082){
cmp=0;
continue;
}
//lecture des donnees SAFRAN, ie raster incomplet
//sprintf(filename,"%s%i",pathData,i);
sprintf(inputDataSafran,"%s/%i/TAIR_%i",PATH_TO_DATA_SAFRAN,numA,step);
FILE *fTempSafran=fopen(inputDataSafran,"rb");
if (!fTempSafran)
continue;
int rfread=fread(pValT, sizeof(double), NB_PT_SAFRAN, fTempSafran);
if (rfread!=NB_PT_SAFRAN){
printf("can't read %i double in file %s, read only %i\n",NB_PT_SAFRAN,inputDataSafran,rfread);
return;
}
fclose(fTempSafran);
cmp++;
for(int locSafran=0;locSafran<NB_PT_SAFRAN;locSafran++){
pValT[locSafran]+=T_ZERO_OBSOLU;
pTday[locSafran]+=pValT[locSafran];
pIjxy[locSafran]+=rhofct(pValT[locSafran]);
}
if (!(step%PERIOD_HEURE)){
for(int locSafran=0;locSafran<NB_PT_SAFRAN;locSafran++){
//if (pTday[locSafran]/24>10)
// printf("%e %e\n",pTday[locSafran]/24,pIjxy[locSafran]);
if (cmp==PERIOD_HEURE){
fprintf(foutput,"%lf %lf\n",pTday[locSafran]/(1.0*PERIOD_HEURE),pIjxy[locSafran]/(1.0*PERIOD_HEURE));
if (pIjxy[locSafran]/(1.0*PERIOD_HEURE)>Vmax)
printf("BUG %i\n",step);
}
pTday[locSafran]=0;
pIjxy[locSafran]=0;
}
cmp=0;
}
}
}
free(pTday);
free(pIjxy);
free(pValT);
fclose(foutput);
}
int main(){
FILE *pf=fopen("rho.txt","w");
for(int ii=-20;ii<45;ii++)
fprintf(pf,"%e %e\n",1.0*ii,rhofct(1.0*ii));
fclose(pf);
safranToIxyTday();
printf("%e %e %e\n",Tmin,Vmin,Vmax);
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment