#include <cstdlib>
#include <iostream>
#include <math.h>
using namespace std;
float geopot[200][200];
float pi,cuento,d,f0,g,X[200][200];
float nf=1;
int files,cols,equador;
void filtre2();
void filtre();
float coriolis (int fila,float dy,float y0,int
doshemisferis );
float derivada2 (float x1,float x2,float x3,float
x4,int direccio,int fila);
float derivada0 (float x1,float x2,float x3,float
x4,float x5,int fila,int col,int direccio);
float derivada (float x1,float x2,int direccio,int
fila);
float dx (int fila);
float valorabsolut(float residu);
char res,enter;
int main(int argc, char *argv[])
{
float
vort[200][200],f,gradient[2][200][200],vgeostrofic[2][200][200],resultatinicial[200][200];
float R[200][200],dp,Fs[200][200],residuacumulat,laplaz,final,advecvort;
int
fila,col,fo,comptador,hora,direccio;
float div,z[100000];
float dades[10];
float
dy,ddx,x0,y0;
int
res,k,j,linia,lcadena,doshemisferis;
char
resposta1,car,cadena[500];
int
sortida1,sortida3,resposta2,pbaixos,fcorrec;
float sortida2;
int horafinal,interval;
FILE *file2 =NULL;
FILE *file1=NULL;
FILE *file3=NULL;
file2=fopen("zeta500.dat","w");
file1=fopen("zeta500.txt","r");
file3= fopen("zeta500.ctl","r+");
g=9.81;
d=55600;
pi=3.14159;
//llegim
la informaciķ de la capcelera del fixer: XDEF,x0,dx,YDEF,y0,dy
for (k=1;k<=6;k++)
fscanf(file1,"%f\n",&dades[k]);
cols=(int)dades[1];
files=(int)dades[4];
ddx=dades[3];
x0=dades[2];
y0=dades[5];
dy=dades[6];
d=d*dy*2;
k=j=linia=0;
//Actualitzem
el fitxer zeta500.ctl amb la informaciķ anterior
do
{
fread(&car,1,1,file3);
k++;
if (linia==3)
j++;
if (car==10)
linia++;
}
while (linia<4);
sprintf(cadena,"%d LINEAR %.1f
%.1f\n",cols,x0,ddx);
fseek(file3,k-j+8,SEEK_SET);
fwrite(&cadena,strlen(cadena),1,file3);
lcadena=strlen(cadena);
sprintf(cadena,"YDEF
%d LINEAR %.1f %.1f\nZDEF 1 LEVELS 500\nTDEF 900 LINEAR 12Z12mar2001 1hr\nVARS
1\nZ 1 99 altura geopotencial de 500 hPa\nENDVARS\n
",files,y0,dy);
fseek(file3,k-j+9+lcadena,SEEK_SET);
//enter=13;
//fwrite(&enter,1,1,file3);
fwrite(&cadena,strlen(cadena),1,file3);
//comencem diāleg per a entrar dades
system("cls");
printf("\ninterval de temps, en minuts ?\n\n\t1 -
15\n\t2 - 30\n\t3 - 60\n\tenter - 15\n");
fflush(stdin);
resposta1=getchar();
switch(resposta1)
{
case 49:
interval=15;
break;
case 50:
interval=30;
break;
case 51:
interval=60;
break;
default:
interval=15;
break;
}
printf("has triat %d minuts ",interval);
resposta1='s';
do
{
printf("\nquants d'intervals vols calcular\? ");
scanf("%d",&horafinal);
if (horafinal>100)
{
printf("\nhas triat %d intervals; es correcte ? (s/n)",horafinal);
fflush(stdin);
resposta1=getchar();
}
printf("has triat %d intervals",horafinal);
}
while (resposta1!='s');
printf("\nvelocitat de rotacio de la terra?\n\n\t1 -
1x\n\t2 - 0.5x\n\t3 - 0.25x\n\t4 - 2x\n\t5 - 4x\n\tenter - 1x\n");
fflush(stdin);
resposta2=(int)getchar();
printf("has
triat una velocitat de rotacio de
");
switch(resposta2)
{
case 49:
nf=1;
printf(" 1x");
break;
case 50:
nf=0.5;
printf(" 0.5x");
break;
case 51:
nf=0.25;
printf(" 0.25x");
break;
case 52:
nf=2;
printf(" 2x");
break;
case 53:
nf=4;
printf(" 4x");
break;
default:
nf=1;
printf(" 1x");
break;
}
fcorrec=0;
do
{
printf("\nles dades son de nomes un hemisferi?
(s/n) ");
fflush (stdin);
resposta1=getchar();
}
while ((resposta1 !='s')&&(resposta1 !='n'));
if (resposta1=='s')
{
doshemisferis=0;
do
{ printf("vols treure el filtre
passa-baixos? (s/n) ");
fflush(stdin);
resposta1=getchar();
}
while ((resposta1
!='s')&&(resposta1 !='n'));
if (resposta1=='s')
pbaixos=0;
else
pbaixos=1;
}
else
{
doshemisferis=1;
do
{ printf("vols posar el filtre
passa-baixos? (s/n) ");
fflush(stdin);
resposta1=getchar();
}
while ((resposta1
!='s')&&(resposta1 !='n'));
if (resposta1=='s')
pbaixos=1;
else
pbaixos=0;
do
{ printf("\nvols canviar el signe del parametre de
coriolis? (s/n)
");
fflush(stdin);
resposta1=getchar();
}
while
((resposta1 !='s')&&(resposta1 !='n'));
if
(resposta1=='s')
fcorrec=1;
else
fcorrec=0;
do
{
printf("\nvols alterar el valor d'f a la
zona equatorial? (s/n)
");
fflush(stdin);
resposta1=getchar();
}
while ((resposta1
!='s')&&(resposta1 !='n'));
if (resposta1=='s')
equador=1;
else
equador=0;
}
//filtrem
les dades del fitxer, ja que de vegades en falten!
fscanf(file1,"%f\n",&z[1]);
for (k=2;k<=(files*cols);k++)
{
fscanf(file1,"%f\n",&z[k]);
if
(fabs((z[k]-z[k-1])>1000))
z[k]=z[k-1];
}
//llegim
la informaciķ del fitxer zeta500.txt
k=1;
for (fila=1;fila<=files;fila++)
{
for(col=1;col<=cols;col++)
{ geopot[fila][col]=z[k];
sortida1=geopot[fila][col];
sortida2=(float)sortida1;
fwrite(&sortida2,4,1,file2);
k++;
}
}
//comenįa
el programa de cālcul de la tendčncia
hora=1;
do
{
switch (doshemisferis)
{ case 1:
if
(pbaixos==0)
{ if
(comptador>3000)
filtre();
else if
(fabs(fila-files/2)<5)
filtre();
}
else
filtre();
break;
case 0:
if
(pbaixos==1)
filtre();
break;
}
for(fila=2;fila<=files-1;fila++)
{ dp=dx(fila);
for(col=2;col<=cols-1;col++)
{ f=coriolis(fila,dy,y0,doshemisferis);
vort[fila][col]=(g/(f*dp*d))*(geopot[fila-1][col]+geopot[fila+1][col]+geopot[fila][col-1]+geopot[fila][col+1]-4*geopot[fila][col])+f
;
}
}
//ja tenim la vorticitat absoluta a la matriu vort
// ara anem a calcular el gradient de la vorticitat
//gradient cap a l'est
for(fila=2;fila<=files-1;fila++)
{
for(col=2;col<=cols-1;col++)
{ direccio=0;
if((col==2)||(col==cols-1)||(fila==2)||(fila==files-1))
gradient[0][fila][col]=derivada0(vort[fila][col-1],vort[fila][col+1],vort[fila-1][col],vort[fila+1][col],vort[fila][col],fila,col,direccio);
else if
((col==3)||(col==cols-2))
gradient[0][fila][col]=derivada(vort[fila][col-1],vort[fila][col+1],direccio,fila);
else
gradient[0][fila][col]=derivada2(vort[fila][col-2],vort[fila][col-1],vort[fila][col+1],vort[fila][col+2],direccio,fila);
}
}
//gradient cap al nord
for(col=2;col<=cols-1;col++)
{ for(fila=2;fila<=files-1;fila++)
{ direccio=1;
if((col==2)||(col==cols-1)||(fila==2)||(fila==files-1))
gradient[1][fila][col]=derivada0(vort[fila][col-1],vort[fila][col+1],vort[fila-1][col],vort[fila+1][col],vort[fila][col],fila,col,direccio);
else
if
((fila==3)||(fila==files-2))
gradient[1][fila][col]=derivada(vort[fila-1][col],vort[fila+1][col],direccio,fila);
else
gradient[1][fila][col]=derivada2(vort[fila-2][col],vort[fila-1][col],vort[fila+1][col],vort[fila+2][col],direccio,fila);
}
}
//ara
anem a calcular el vent geostrōfic
//component
cap al nord del vent geostrōfic
for
(fila=2;fila<=files-1;fila++)
{ f=coriolis(fila,dy,y0,doshemisferis);
for(col=2;col<=cols-1;col++)
{ direccio=0;
if((col==2)||(col==cols-1)||(fila==2)||(fila==files-1))
vgeostrofic[1][fila][col]=(g/f)*derivada0(geopot[fila][col-1],geopot[fila][col+1],geopot[fila-1][col],geopot[fila+1][col],geopot[fila][col],fila,col,direccio);
else
{ if((col==3)||(col==cols-2))
vgeostrofic[1][fila][col]=(g/f)*derivada(geopot[fila][col-1],geopot[fila][col+1],direccio,fila);
else
vgeostrofic[1][fila][col]=(g/f)*derivada2(geopot[fila][col-2],geopot[fila][col-1],geopot[fila][col+1],geopot[fila][col+2],direccio,fila);
}
}
}
//component cap a l'est del vent geostrōfic
for (col=2;col<=cols-1;col++)
{ for(fila=2;fila<=files-1;fila++)
{
f=coriolis(fila,dy,y0,doshemisferis);
direccio=1;
if((col==2)||(col==cols-1)||(fila==2)||(fila==files-1))
vgeostrofic[0][fila][col]=-(g/f)*derivada0(geopot[fila][col-1],geopot[fila][col+1],geopot[fila-1][col],geopot[fila+1][col],geopot[fila][col],fila,col,direccio);
else
if((fila==3)||(fila==files-2))
vgeostrofic[0][fila][col]=-(g/f)*derivada(geopot[fila-1][col],geopot[fila+1][col],direccio,fila);
else
vgeostrofic[0][fila][col]=-(g/f)*derivada2(geopot[fila-2][col],geopot[fila-1][col],geopot[fila+1][col],geopot[fila+2][col],direccio,fila);
}
}
//ara
anem a fer el producte escalar del v. geostrōfic pel gradient de la vorticitat
for(col=2;col<=cols-1;col++)
{ for
(fila=2;fila<=files-1;fila++)
{
f=coriolis(fila,dy,y0,doshemisferis);
advecvort=-((vgeostrofic[0][fila][col])*gradient[0][fila][col]+(vgeostrofic[1][fila][col])*gradient[1][fila][col]);
if
((advecvort>0)&&(nf==1)&&(doshemisferis==1))
advecvort=advecvort*1.05; // potenciem l'advecciķ de vorticitat +
switch (fcorrec)
{
case 1:
if
(fila<((int)(files/2)))
X[fila][col]=-(f/g)*advecvort;
else
X[fila][col]=(f/g)*advecvort;
break;
default:
X[fila][col]=(f/g)*advecvort;
}
}
}
//ara
anem a fer l'algoritme de la relaxaciķ simultānia.
comptador=0;
for(col=2;col<=cols-1;col++)
for(fila=2;fila<=files-1;fila++)
Fs[fila][col]=X[fila][col];
// ja tenim la laplaciana;
//llimem la frontera
for (fila=1;fila<=files;fila++)
{
X[fila][1]=X[fila][2];
X[fila][cols]=X[fila][cols-1];
}
for(col=1;col<=cols;col++)
{
X[1][col]=X[2][col];
X[files][col]=X[files-1][col];
}
do
{
residuacumulat=0;
for(col=2;col<=cols-1;col++)
{
for(fila=2;fila<=files-1;fila++)
{ dp=dx(fila);
R[fila][col]=X[fila+1][col]+X[fila-1][col]+X[fila][col+1]+X[fila][col-1]-4*X[fila][col]-d*d*Fs[fila][col];
residuacumulat=residuacumulat+fabs(R[fila][col]);
}
}
for(col=2;col<=cols-1;col++)
for(fila=2;fila<=files-1;fila++)
X[fila][col]=X[fila][col]+0.25*R[fila][col];
comptador++;
if
(comptador>20000)
{ printf("El proces no convergeix\n\n");
printf("files=%d,
columnes=%d\n",files,cols);
printf("El
residu acumulat es : %f\n",residuacumulat);
fclose(file1);
fclose(file2);
system("PAUSE");
return
EXIT_SUCCESS;
}
}
while(residuacumulat>0.001);
system("cls");
printf("Comptador=%d\nInterval=%d\n\n",comptador,hora);
switch (doshemisferis)
{ case 1:
if
(pbaixos==0)
{ if
(comptador>3000)
filtre2();
else if
(fabs(fila-files/2)<5)
filtre2();
}
else
filtre2();
break;
case 0:
if
(pbaixos==1)
filtre2();
break;
}
//llimem la frontera
for (fila=1;fila<=files;fila++)
{
X[fila][1]=X[fila][2];
X[fila][cols]=X[fila][cols-1];
}
for(col=1;col<=cols;col++)
{
X[1][col]=X[2][col];
X[files][col]=X[files-1][col];
}
for(fila=1;fila<=files;fila++)
for(col=1;col<=cols;col++)
{ switch
(doshemisferis)
{ case
1:
if
(fabs(fila-files/2)>0)
geopot[fila][col]=geopot[fila][col]+60*interval*X[fila][col]-60*interval*0.0006;
break;
case 0:
geopot[fila][col]=geopot[fila][col]+(60*interval)*(X[fila][col]);
break;
}
sortida1=geopot[fila][col];
sortida2=(float)sortida1;
fwrite(&sortida2,4,1,file2);
}
hora++;
}
while
(hora<horafinal);
fclose(file1);
fclose(file2);
fclose(file3);
system("PAUSE");
return
EXIT_SUCCESS;
}
void filtre2()
{
int fila,col;
for(col=1;col<=cols;col++)
for
(fila=2;fila<=files-1;fila++)
X[fila][col]=(X[fila-1][col]+2*X[fila][col]+X[fila+1][col])/4;
for (fila=1;fila<=files;fila++)
for(col=2;col<=cols-1;col++)
X[fila][col]=(X[fila][col-1]+2*X[fila][col]+X[fila][col+1])/4;
}
void filtre()
{
int fila,col;
for(col=1;col<=cols;col++)
for
(fila=2;fila<=files-1;fila++)
geopot[fila][col]=(geopot[fila-1][col]+2*geopot[fila][col]+geopot[fila+1][col])/4;
for (fila=1;fila<=files;fila++)
for(col=2;col<=cols-1;col++)
geopot[fila][col]=(geopot[fila][col-1]+2*geopot[fila][col]+geopot[fila][col+1])/4;
}
float coriolis(int fila,float dy,float y0,int
doshemisferis)
{
float lat,f,y;
//lat=(24.5+0.5*fila)*pi/180;
lat=((y0-dy)+dy*fila)*pi/180;
switch (doshemisferis)
{ case
1:
if
(equador==1)
{ if
(((fila-files/2)<3)&&((fila-files/2)>=0))
lat=((y0-dy)+dy*(fila+3))*pi/180;
else if (((fila-files/2)>-3)&&((fila-files/2)<0))
lat=((y0-dy)+dy*(fila-3))*pi/180;}
break;
default:
break;
}
f=2*sin(lat)*2*pi/(24*36);
return(nf*f/100);
}
float derivada2(float x1,float x2,float x3,float
x4,int direccio,int fila)
{
float dp,dist;
dp=dx(fila);
dist = d;
if (direccio==0)
dist=0;
dist=d;
return(x1-8*x2+8*x3-x4)/(12*dist);
}
float derivada0(float x1,float x2,float x3,float
x4,float x5,int fila,int col,int direccio)
{
float a,b,dp,dist;
dp=dx(fila);
dist=d;
if (direccio==1)
{ switch
(fila)
{ case 2:
a=x4;
b=x5;
break;
case 92:
a=x5;
b=x3;
break;
}
}
else
{ dist=dp;
switch(col)
{ case 2:
a=x2;
b=x5;
break;
case 132:
a=x5;
b=x1;
break;
}
}
dist = d;
return((a-b)/dist);
}
float derivada(float x1,float x2,int direccio,int
fila)
{
float dp,dist;
dp=dx(fila);
if
(direccio==0)
dist=dp;
dist=d;
return((x2-x1)/(2*dist));
}
float dx ( int
fila)
{
float lat,dist;
return(d);
}
float valorabsolut(float residu)
{
float valor;
if (residu<0)
valor=-residu;
else
valor=residu;
return(valor);
}