cycle.cpp

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

//************************************************************************************************
//***************  本程序用来计算平流层的气球的航线，目的通过控制气球的高度，使气球在
//***************  一定垂直区域内实现循环运动，使用的资料是17层的间隔6小时的t106资料。
//************************************************************************************************
//-------------------------------------------------------------
//------航线规划主程序开始。
//-------------------------------------------------------------
//--适用于标准层次(10,20,30,50,70,100,150,200,250,300,400,500,700,850hpa)飞行。
//    implicit none
include "cycle.h"

const int	ntime=9,nlevel=17,j7=27,i7=37,t7=49,varnum=6,k7=153;
const float rd=2.871;	   // 定义时间和层次变量
const float	nu=1.49e-5;				//
const float pi=3.14159;				//
int			nz7(i7,j7,nlevel,ntime),nu7(i7,j7,nlevel,ntime),nv7(i7,j7,nlevel,ntime),
			nt7(i7,j7,nlevel,ntime);	        // 定义数组，存储y原始t106资料
float		nz24(i7,j7,nlevel,t7),nu24(i7,j7,nlevel,t7),nv24(i7,j7,nlevel,t7),
			nt24(i7,j7,nlevel,t7);           // 定义数组，存储插值后的t106逐时资料
double		nz24(i7,j7,nlevel,ntime),nu24(i7,j7,nlevel,ntime),nv24(i7,j7,nlevel,ntime),
			nt24(i7,j7,nlevel,ntime);	    // 定义数组，存储t106资料
float		blnradius,uf,vf;
float		profile[varnum,nlevel];



void balloon_17()
{
	int i,l1,l2,l3,l4,l5;									//	 定义整型变量
	int nflyhour,maxhour,nprnttime,minprnt,lbrsn_time;		// 定义飞行时间，最大飞行时间等变量 
	float dlon,dlat,blnrads		 //	 定义经纬度、气球半径变量
	char t106dir[120],date[120],stime[120],profilename[120],myblnroute[120]	// 定义字符型变量，t106路径、日期等
	float verspd[16],distance      //	 垂直速度数组
//  读取设置信息
	FILE *m_file;
	CString m_workfile = "blnrouter.ini";
	if( (m_file  = fopen( m_workfile, "r" )) != null )
	{
		fscanf(m_file,"%s\n",&t106dir);
		fscanf(m_file,"%s\n",&date);
		fscanf(m_file,"%s\n",&stime);
		fscanf(m_file,"%f\n",&dlon);
		if(dlon<60.0 || dlon>150.0)
		{
			AfxMessageBox("经度超出东经60°～150°的范围");
			return;
		}
		fscanf(m_file,"%f\n",&dlat);
		if(dlat<15.0 || dlat>80.0)
		{
			AfxMessageBox("纬度超出北纬15°n～80°n的范围");
			return;
		}
		fscanf(m_file,"%d\n",&nflyhour);
		fscanf(m_file,"%d\n",&lbrsn_time);
		fscanf(m_file,"%s\n",&profilename);
		fscanf(m_file,"%s\n",&myblnroute);
		fscanf(m_file,"%d\n",&nprnttime);
		fscanf(m_file,"%f\n",&blnrads);
		fscanf(m_file,"%f\n",&distance);
		fclose( m_file );
	}
	m_workfile = "vertical_speed.ini";
	if( (m_file  = fopen( m_workfile, "r" )) != null )
	{
		for(int i=1;i<=16;i++)
			fscanf(m_file,"%f\n",&verspd[i]);
		fclose( m_file );
	}

//  求取字符串长度
	l1=strlen(t106dir);				  
	l2=strlen(date);				  
	l3=strlen(profilename);			  
	l4=strlen(myblnroute);			  
	l5=strlen(stime);				  

	maxhour=nflyhour;
	minprnt=nprnttime;
// 调用读取t106原始资料（间隔6小时）的子程序  
	read_t106_china(t106dir(1:l1);
// 调用读取将t106原始资料插值成逐时资料的子程序
	readt10624();
// 调用计算释放站点0～30km（1000hpa～10hpa）廓线资料的子程序
	profile_17(dlon,dlat,lbrsn_time,profilename(1:l3));
// 调用计算气球航线的子程序的子程序
	bln_trace(dlon,dlat,maxhour,minprnt,lbrsn_time,myblnroute(1:l4),blnrads,verspd,distance);
}
//------------------  航线规划主程序结束。 -------------------


//--------------计算气球航线的子程序-------------------------
void bln_trace(float dlon,float dlat,int maxhour,int minprnt,int lbrsn_time,char * bln_router_file_name,
	                    float  blnrads,float verspd[nlevel-1],float distance)
{
	int ht,k,nt,mt;	                // 定义循环变量等
	int nlon,nlat,istep; 				// 定义变量
	int ifprnt;					// 定义变量
	const float dt = 0.1,t0 = 0.0,x0 = 0.0,y0 = 0.0;
	float z0,t,x,y,z;
// 资料范围为：60e～150e（格点24～60），15n～80n（格点6～32）
    float rr;
	float ylat,xlon	    //
	float u,v,w,dz,fx,fy		    // uf,vf为环境风场的分量
//------------------------------------------------------------------------------------------------
// 获取输出文件名。
	if( (m_file  = fopen( bln_router_file_name, "r" )) != null )
	{
		blnradius=blnrads;		//气球半径。
		// 计算释放点所在格点位置。
		nlon=int((dlon-60.)/2.5)+1;		    // 减去23，为了和数组的下标i保持一致
		nlat=int((80.-dlat)/2.5)+1;		    // 为了和数组的下标j保持一致
		// 赋初值
		z0=0.0;	u =0.0;	v =0.0;	uf=0.0;	vf=0.0;
		// 赋时间和局地坐标初始值。
		t=t0;		 //	积分时间
		x=x0;		 //	飞行纬向距离
		y=y0;		 //	飞行经向距离
		z=z0;		 //	上升高度
		istep=0; 	 //  时步数统计变量初值。
		nt=0;         //  小时数统计变量初值。
		rr=6371*1000; //  地球半径：单位为米。
		ifprnt=int(minprnt*60/dt);  //打印间隔。
		istep=0;
		//--------飞行过程-----------------------------------------------------------------------
		while( sqrt(x*x+y*y) <= distance) // 总飞行时间控制积分时间,飞行400km就完成
		{
			ht=nt +lbrsn_time;			          // 时次
			ylat=y/rr*180./pi;  					  //	飞行的纬度距离
			//  判断范围
			if((dlat+ylat)<15.0 || (dlat+ylat)>80.0)
				return;
			xlon=x/(rr*cos(dlat*pi/180.))*180./pi; //	飞行的经度距离
			// 判断范围
			if((dlon+xlon)<60.0 || (dlon+xlon)>150.0)
				return;
			// 由球面坐标计算格点坐标。
			nlon=int((dlon+xlon-60.)/2.5)+1;
			nlat=int(80.-dlat-ylat)/2.5+1;
			// 判断气球飞行的高度	
			if (z>=profile[2,17])
			{
				// 预定高度处为平飞过程 
				uf=  nu24(nlon,nlat,17,ht)/10.; // 预定高度处x方向风速
				vf=  nv24(nlon,nlat,17,ht)/10.; // 预定高度处y方向风速
				w = 0.0;
			}
			// 到达预定高度，垂直速度为0，不再上升
			else if(z < profile[2,17])
			{
				// 预定高度下为上升过程
				if(z<profile[2,1])
				{
					uf=  nu24(nlon,nlat,1,ht)/10.; // 预定高度以下x方向风速
					vf=  nv24(nlon,nlat,1,ht)/10.; // 预定高度以下y方向风速
					w = verspd(1);				  // 预定高度以下的垂直升速
				}
				else
				{
					for(k=1;k<nlevel-1;k++)
					{
						if (z>profile[2,k] && z<profile[2,k+1])
						{
							// 为初始风速风向赋值，进行插值并转换成实际值。
							dz=1.0*(z-nz24(nlon,nlat,k,ht))/(nz24(nlon,nlat,k+1,ht)-nz24(nlon,nlat,k,ht));
							uf=(dz*(nu24(nlon,nlat,k+1,ht)-nu24(nlon,nlat,k,ht))+nu24(nlon,nlat,k,ht))/10.;
							vf=(dz*(nv24(nlon,nlat,k+1,ht)-nv24(nlon,nlat,k,ht))+nv24(nlon,nlat,k,ht))/10.;
							w = verspd(k);
						}
					}
				}
			}
			// 调用kutta子程序积分。
			kutta(x,y,z,u,v,w,t,dt,fx,fy);
			// 计算小时数
			nt=int(istep*dt/3600.);
			// 计算分钟数
			mt=int(istep*dt/60.);
			// 到达时间间隔输出地理坐标值一次。
			if((istep/10)*10 == 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;
			// 准备结束程序。
		}
		
		//--------- 距离释放点到达预定范围时，开始下降到西风带200hpa，准备返回
		while(z>=profile[2,10])
		{
			w = -4.;
			ht = nt +lbrsn_time;			          // 时次
			ylat = y/rr*180./pi;  					  //	飞行的纬度距离
			xlon = x/(rr*cos(dlat*pi/180.))*180./pi; //	飞行的经度距离
			// 由球面坐标计算格点坐标。
			nlon = int((dlon+xlon-60.)/2.5)+1;
			nlat = int(80.-dlat-ylat)/2.5+1;
			// 判断气球飞行的高度	
			if (z >= profile[2,17])
			{
				// 预定高度处为平飞过程 
				uf=  nu24(nlon,nlat,17,ht)/10.; // 预定高度处x方向风速
				vf=  nv24(nlon,nlat,17,ht)/10.; // 预定高度处y方向风速
			}
			else
			{
				for(k=9;k<=nlevel-1;k++)
				{
					if (z>profile[2,k] && z<profile[2,k+1])
					{
						// 为初始风速风向赋值，进行插值并转换成实际值。
						dz=1.0*(z-nz24(nlon,nlat,k,ht))/(nz24(nlon,nlat,k+1,ht)-nz24(nlon,nlat,k,ht));
						uf=(dz*(nu24(nlon,nlat,k+1,ht)-nu24(nlon,nlat,k,ht))+nu24(nlon,nlat,k,ht))/10.;
						vf=(dz*(nv24(nlon,nlat,k+1,ht)-nv24(nlon,nlat,k,ht))+nv24(nlon,nlat,k,ht))/10.;
					}
				}
			}
			// 调用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;
		}
		//--------------------------------------------------------------------------------------
		//------ 下降到预定高度200hpa后，停止下降，开始返回到出发点附近（同一经度上）------ 
		//--------------------------------------------------------------------------------------
		while(((dlon+xlon)-dlon) < 0.05)
		{
			ht=nt +lbrsn_time;			          // 时次
			ylat=y/rr*180./pi;  					  //	飞行的纬度距离
			xlon=x/(rr*cos(dlat*pi/180.))*180./pi; //	飞行的经度距离
			
			// 由球面坐标计算格点坐标。
			nlon=int((dlon+xlon-60.)/2.5)+1;
			nlat=int(80.-dlat-ylat)/2.5+1;
			// 预定高度处为平飞过程 
			uf =  nu24(nlon,nlat,10,ht)/10.; // 预定高度处x方向风速
			vf =  nv24(nlon,nlat,10,ht)/10.; // 预定高度处y方向风速