cycle.cpp

气象数据计算超等许,能够计算在风厂中气象数据的票性诡计

			w = 0.0;                        // 到达预定高度，垂直速度为0，不再上升
			// 调用kutta子程序积分。
			kutta(x,y,z,u,v,w,t,dt,fx,fy);
			// 计算小时数
			nt=int(istep*dt/3600.);
			// 计算分钟数
			mt=int(istep*dt/60.);
			// 到达时间间隔输出地理坐标值一次。
			if((istep/ifprnt)*ifprnt == istep)
			{
				fprintf(m_file,"%d\n",nt);
				fprintf(m_file,"%d\n",mt+1);
				fprintf(m_file,"%f\n",dlon+xlon);
				fprintf(m_file,"%f\n",dlat+ylat);
				fprintf(m_file,"%f\n",z);
				fprintf(m_file,"%f\n",u);
				fprintf(m_file,"%f\n",v);
				fprintf(m_file,"%f\n",uf);
				fprintf(m_file,"%f\n",vf);
				fprintf(m_file,"%f\n",w);
			}
			// 积分步数累计，置于此处是为了当istep＝0时，也输出一次结果				   
			istep=istep+1;
			//判断是否回到原来的出发点或附近
		}
		fclose( m_file );
	}
}
///////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////
//函数fx(x,y,u,v,t)用于定义x方向的力。
float fx(float x,float y,float u,float v,float t)  
{
	float fx,wr,cd;
	wr=sqrt((uf-u)*(uf-u)+(vf-v)*(vf-v));
	if (blnradius <= 0.0)
		fx=0.0;
	else
		fx=1.5/blnradius*cd(wr)*(uf-u)*wr;
	return fx;
}
//-------------------------------------------------------------
// 函数fy(x,y,u,v,t)用于定义y方向的力。
//-------------------------------------------------------------
float fy(float x,float y,float u,float v,float t)
{
	float fy,wr,cd;
	wr=sqrt((uf-u)*(uf-u)+(vf-v)*(vf-v));
	if (blnradius <= 0.0)
		fx=0.0;
	else
		fy=1.5/blnradius*cd(wr)*(vf-v)*wr;
	return fy;
}
//-------------------------------------------------------------
//函数cd(v)用于定义球体的阻力系数。
//-------------------------------------------------------------
float cd(float v)
{
	float re,cd;
	re=abs(v)*(2*blnradius)/nu;
	if(re==0.)						cd=0;
	if(re>0. && re <= 1.)			cd=24./re;
	if(re>1. && re <=  400.)		cd=24./re**0.646;
	if(re>400. && re <=  3.0e+5)	cd=0.5;
	if(re> 3.0e+5 && re <=  2.0e+6)	cd=3.66e-4*re**0.4275;
	if(re > 2.0e+6)					cd=0.18;
	return cd;
}
//-----------定义子例行程序，龙格－库塔四阶积分-----------------
//-------------------------------------------------------
//四阶runger-kutta方法，求解如下形式的二阶常微分议程组。
//h表示时间t在积分时的时间步长。
//f1，f2在其他函数中定义。
//x,y,u,v,t,h,f1,f2均在主程序中说明。
void kutta(float x,float y,float z,float u,float v,float w,float t,float h,float f1,float f2)
{
	float d1x,d1y,d1u,d1v,d2x,d2y,d2u,d2v,d3x,d3y,d3u,d3v,d4x,d4y,d4u,d4v;

	d1x=h*u;
	d1y=h*v;
	d1u=h*f1(x,y,u,v,t);
	d1v=h*f2(x,y,u,v,t);

	d2x=h*(u+d1u/2.);
	d2y=h*(v+d1v/2.);
	d2u=h*f1(x+d1x/2.,y+d1y/2.,u+d1u/2.,v+d1v/2.,t+h);
	d2v=h*f2(x+d1x/2.,y+d1y/2.,u+d1u/2.,v+d1v/2.,t+h);

	d3x=h*(u+d2u/2.);
	d3y=h*(v+d2v/2.);
	d3u=h*f1(x+d2x/2.,y+d2y/2.,u+d2u/2.,v+d2v/2.,t+h);
	d3v=h*f2(x+d2x/2.,y+d2y/2.,u+d2u/2.,v+d2v/2.,t+h);

	d4x=h*(u+d3u);
	d4y=h*(v+d3v);
	d4u=h*f1(x+d3x,y+d3y,u+d3u,v+d3v,t+h);
	d4v=h*f2(x+d3x,y+d3y,u+d3u,v+d3v,t+h);

	t=t+h;
	x=x+(d1x+2.*d2x+2.*d3x+d4x)/6.;
	y=y+(d1y+2.*d2y+2.*d3y+d4y)/6.;
	u=u+(d1u+2.*d2u+2.*d3u+d4u)/6.;
	v=v+(d1v+2.*d2v+2.*d3v+d4v)/6.;

	z=z+w*h;
}
//\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
//\\\\\	本子程序用来计算两天内的逐时的风场u，v；位势z和温度t的值
/////////////////////////////////////////////////////////////////////////////////////////////////
void readt10624()
{
	int i,j,k,nt,ht;	// 定义变量
	int ntime,nlevel; 		    // 定义时间和层次变量
									//资料范围为：60e～150e（格点24～60），15n～80n（格点6～32）
	float ndltz,ndltu,ndltv,ndltt;	        // 定义变化量的变量
	//－计算1-48小时内的逐时风场
	for(nt=1;nt<=ntime-1;nt++)
	{
		for(k=1;k<=nlevel;k++)
		{
			for(i=1;i<=37;i++)
			{
				for(j=1;j<=27;j++)
				{
					//-计算变化量
					ndltz=(nz7[i,j,k,nt+1]-nz7[i,j,k,nt])/6.;
					ndltu=(nu7[i,j,k,nt+1]-nu7[i,j,k,nt])/6.;
					ndltv=(nv7[i,j,k,nt+1]-nv7[i,j,k,nt])/6.;
					ndltt=(nt7[i,j,k,nt+1]-nt7[i,j,k,nt])/6.;
					//-风场时间插值
					for(ht=1;ht<=7;ht++)
					{
						nz24[i,j,k,ht+6*(nt-1)]= nz7[i,j,k,nt]+(ht-1)*ndltz;
						nu24[i,j,k,ht+6*(nt-1])= nu7[i,j,k,nt]+(ht-1)*ndltu;
						nv24[i,j,k,ht+6*(nt-1)]= nv7[i,j,k,nt]+(ht-1)*ndltv;
						nt24[i,j,k,ht+6*(nt-1])= nt7[i,j,k,nt]+(ht-1)*ndltt;
					}
				}
			}
		}
	}
}
//------------------------------------------------------
//--------readt10624子程序结束。
//------------------------------------------------------
////////////////////////////////////////////////////////////////////////////////////////
//-------------------本子程序用来读取17层的t106文件资料----------------------------
void read_t106_china(CString t106file)
{
	int i,j,k,nt,j7,i7,k7;	// 定义变量
	int ntime,nlevel; 		// 定义时间和层次变量
// 获取资料路径
	FILE *m_file;
	if( (m_file  = fopen( t106file, "r" )) != null )
	{
// 读取	
		for(nt=1;nt<=ntime;nt++)
		{
			for(k=1;k<=nlevel;k++)
			{
				for(i=1;i<=i7;i++)
					for(j=1;j<=j7;j++)
						fscanf(m_file,"%s\n",&nz7[i,j,k,nt]);		//  位势高度
				for(i=1;i<=i7;i++)
					for(j=1;j<=j7;j++)
						fscanf(m_file,"%s\n",&nu7[i,j,k,nt]);		//  水平风场u分量
				for(i=1;i<=i7;i++)
					for(j=1;j<=j7;j++)
						fscanf(m_file,"%s\n",&nv7[i,j,k,nt]);		//  水平风场v分量
				for(i=1;i<=i7;i++)
					for(j=1;j<=j7;j++)
						fscanf(m_file,"%s\n",&nt7[i,j,k,nt]); 	    //  温度场（绝对温度）
			}
		}
		fclose( m_file );
	}
}

//\\本子程序用来计算所选站点上空1000hpa～10hpa共17层标准等压面的廓线数据
//////////////////////////////////////////////////////////////////////////////////
void profile_17(float dlon,float dlat,int lbrsn_time,CString t106file)	
{
	int i,k,varnum,j7,i7;				// 定义变量
	int ntime,nlevel; 					// 定义时间和层次变量
	int  nlon,nlat,ppp;					// 定义经纬度格点变量
	//资料范围为：60e～150e（格点24～60），15n～80n（格点6～32）
	// 计算所选站点在数组中的位置
	nlon=int((dlon-60.)/2.5)+1;		    // 减去23，为了和数组的下标i保持一致
	nlat=int((80.-dlat)/2.5)+1;		    // 为了和数组的下标j保持一致
	// 计算廓线资料，按气压p，位势高度，u风场，v风场，温度t和密度顺序存放
	FILE* m_file;
	if( (m_file  = fopen( t106file, "w+r" )) != null )
	{
		//  计算开始释放气球时刻的廓线资料
		i=lbrsn_time+1;
		//  计算各层气压
		for(int k=1;k<nlevel+1;k++)  //
		{
			if(k = 1)
				ppp=1000;
			else if(k = 2)
				ppp=925;
			else if(k = 3)
				ppp=850;
			else if(k = 4)
				ppp=700;
			else if(k = 5)
				ppp=600;
			else if(k = 6)
				ppp=500;
			else if(k = 7)
				ppp=400;
			else if(k = 8)
				ppp=300;
			else if(k = 9)
				ppp=250;
			else if(k = 10)
				ppp=200;
			else if(k = 11)
				ppp=150;
			else if(k = 12)
				ppp=100;
			else if(k = 13)
				ppp=70;
			else if(k = 14)
				ppp=50;
			else if(k = 15)
				ppp=30;
			else if(k = 16)
				ppp=20;
			else if(k = 17)
				ppp=10;
			else
				AfxMessageBox("k的值不在 1～17 之间");
			//  赋值给廓线数组
			profile[1,k]=ppp;				              //	 气压
			profile[2,k]=nz24(nlon,nlat,k,i);	          //	 位势高度
			profile[3,k]=nu24(nlon,nlat,k,i)/10.;	      //	 u风场
			profile[4,k]=nv24(nlon,nlat,k,i)/10.;	      //	 v风场
			profile[5,k]=nt24(nlon,nlat,k,i)/10.-273.15;	  //	 温度（摄氏温度）
			// 计算密度，根据p=nu*r*i计算
			nu=ppp/(rd*nt24(nlon,nlat,k,i)/10.);  
			profile[6,k]=nu; 							  //  密度
		}
		//  写文件题头
		CString str = "气压(hpa)\t位势高度(m)\tu风场(m/s)\tv风场(m/s)\t温度(℃)\t密度(kg/m3)\n";
		fprintf(m_file,"%s\n",str);
		//  写入文件myprofile.dat，输出
		for(k=1;k<=nlevel;k++) 
		{
			fprintf(m_file,"%f\t",profile[1,k]);
			fprintf(m_file,"%f\t",profile[2,k]);
			fprintf(m_file,"%f\t",profile[3,k]);
			fprintf(m_file,"%f\t",profile[4,k]);
			fprintf(m_file,"%f\t",profile[5,k]);
			fprintf(m_file,"%f\n",profile[6,k]);
		}
		AfxMessageBox("廓线计算完成");
		fclose(m_file);
	}
}