ImageVerifierCode 换一换
格式:DOCX , 页数:19 ,大小:337.25KB ,
资源ID:9242394      下载积分:3 金币
快捷下载
登录下载
邮箱/手机:
温馨提示:
快捷下载时,用户名和密码都是您填写的邮箱或者手机号,方便查询和重复下载(系统自动生成)。 如填写123,账号就是123,密码也是123。
特别说明:
请自助下载,系统不会自动发送文件的哦; 如果您已付费,想二次下载,请登录后访问:我的下载记录
支付方式: 支付宝    微信支付   
验证码:   换一换

加入VIP,免费下载
 

温馨提示:由于个人手机设置不同,如果发现不能下载,请复制以下地址【https://www.bingdoc.com/d-9242394.html】到电脑端继续下载(重复下载不扣费)。

已注册用户请登录:
账号:
密码:
验证码:   换一换
  忘记密码?
三方登录: 微信登录   QQ登录  

下载须知

1: 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。
2: 试题试卷类文档,如果标题没有明确说明有答案则都视为没有答案,请知晓。
3: 文件的所有权益归上传用户所有。
4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
5. 本站仅提供交流平台,并不能对任何下载内容负责。
6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

版权提示 | 免责声明

本文(A地区边界层结构特征的数值模拟研究.docx)为本站会员(b****0)主动上传,冰点文库仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对上载内容本身不做任何修改或编辑。 若此文所含内容侵犯了您的版权或隐私,请立即通知冰点文库(发送邮件至service@bingdoc.com或直接QQ联系客服),我们立即给予删除!

A地区边界层结构特征的数值模拟研究.docx

1、A地区边界层结构特征的数值模拟研究A地区边界层结构特征的数值模拟研究文章利用2019年A气象站的常规观测资料和探空资料,对A盆地地区的边界层结构特征进行系统的分析,并利用WRF模式的两种边界层参数化方案(YSU方案和MYJ方案),分别选取2019年的7月56日和1月56日作为夏季和冬季的个例,对A盆地地区的大气边界层结构特征进行了数值模拟。结果表明:无论夏季和冬季,WRF模式在A地区均有较好的适用性;夏冬两季的7:00只存在不稳定和中性两种状态,而19:00则只存在稳定和中性这两种状态;夏季的边界层高度普遍高于冬季,19:00的高度高于7:00,WRF模式的YSU方案能够较好的模拟边界层的高度

2、;对于位温和比湿,通过分析,可以具体将其廓线类型分为四种类型,WRF对夏季位温廓线模拟效果较好,而在冬季,对比湿廓线的模拟效果比较理想,这可能与冬季湿度小,天气稳定有关;感(潜)热通量在白天均为正值,晚上为负值,夏季地面向大气传输能量的强度和时间明显强于冬季,在两个季节中,WRF模式模拟的的感热通量的变化趋势与实测值均基本吻合,但是对于潜热通量,冬季的模拟效果明显优于夏季。关键词:A盆地;边界层结构特征;WRF模式;边界层参数化方案;适用性第一章 引言1.1 研究背景及意义大气边界层,是指受地球表面摩擦以及热过程和蒸发显著影响的大气层,即大气底部受地面影响最大的一层【1-2】。大气边界层的平均

3、厚度约为地表以上1000m范围【3】,地表提供的物质和能量主要消耗和扩散在大气边界层内。边界层大气受地球表面的影响,对表面强迫力的响应时间约为1小时或更短,因此大气边界层与人类活动的关系最密切、最直接。大气污染物的散布与其密切相关,大气边界层内的大气扩散能力和其稳定度密切相关,大气稳定度随着气温层结的分布而变化,是直接影响大气污染物扩散的重要因素,大气边界层的高度是空气污染预报的重要指标之一【4-6】。另外,大气边界层是地球大气动量、热量和各种物质上下输送的重要通道,因此,大气边界层中发生的物理过程是大尺度天气气候形成的重要因子,也是天气气候影响人类生活以及人类活动影响大气的通道【7】。大气边

4、界层的物理过程与下垫面的物理性质具有紧密关系,由于边界层过程又受到天气尺度背景场的制约与影响,并且也对天气尺度场起着反馈作用,因此两者关系十分密切。大气边界层是整个地球大气层中最复杂,也是最关键的部分之一,大气边界层起着大气和地球之间的桥梁作用,它们之间的物质和能量交换都是通过大气边界层过程来实现【8】。为提高天气气候预报的准确性,预测大气污染以及推动其他应用气象学科与气候学的发展,深入研究大气边界层的结构特征是非常有必要的。在我国西北半干旱地区,主要为黄土高原农牧交错区,这里既是气候变化敏感区,又是生态环境脆弱带,还是黄河中上游水土保持重点区域。该区域是雨养农业区,农牧林业生产和生态环境对气

5、候条件的依赖性极强,气候变化对社会经济影响重大【9】。近些年,很多学者对我国西北干旱区的边界层特征进行了研究【1013】,相比之下,对半干旱区的边界层特征的研究相对较少,而且存在地区区域性不强、特征参数不够具体等一系列问题,缺乏系统的分析。随着数值模式的推广,WRF模式已经成为科研和业务中的一项重要技术。目前,国内应用WRF模式主要是用来研究典型天气过程【1417】以及对城市边界层结构的模拟【1820】,而对于WRF模式在半干旱地区的应用问题的研究相对较少。文章将系统的分析A地区的边界层结构特征,并使用WRF模式对处于半干旱区的A盆地地区的夏冬两季的边界层结构进行模拟,通过模拟值与观测资料的对

6、比,探究WRF模式在该地区的适用性,进而为半干旱地区边界层结构特征的研究提供一种参考。1.2资料选取文章所用的资料为:2019年A气象站的探空资料、常规观测资料,2019年兰州大学半干旱气候和环境观测站(SACOL)站的通量观测数据以及2019年NCEP11再分析资料。1.3 方法介绍1.1.1 Turner稳定度分类方法 我们将表征大气稳定性质的度量称为大气稳定度,其可以分为稳定性、中性和不稳定性三种类型。所谓稳定性,是指周围大气温度和湿度的铅直分布,具有使受扰气块回到原来位置的趋势;若使受扰气块既无回到原来位置又无远离原来位置的趋势,则称大气为中性的;若使受扰气块具有继续远离原来位置的趋势

7、,则称大气具有不稳定性。对于大气稳定度的分类,有很多方法,本文采用的是Turner稳定度分类法(简称T方法),将大气稳定度具体分为A、B、C、D、DE、E和F七个等级(A代表极不稳定,B代表中等不稳定,C代表弱不稳定、D代表中性、DE类的稳定性介于中性和弱稳定之间,定义为微弱稳定,E代表弱稳定,F代表中等稳定),其中A、B、C属于不稳定性,D属于中性,DE、E、F属于稳定性。T方法需要使用地面常规观测资料中的云高、总云量、低云量和10m高风速并利用公式 (1)其中为地理纬度,为赤纬,为时角,计算出太阳高度角h,进而再结合日射等级确定规则表以及Turner稳定度分级方法表即可对稳定度进行分级。1

8、.3.2综合判定法a位温廓线法根据探空资料中的温度资料以及气压资料,利用公式 (2)计算位温数据。在位温廓线图上,通常将边界层内位温梯度明显不连续的高度作为边界层高度。如在白天对流发展时,就将对流混合层顶(上部稳定层的低部)作为边界层高度;而夜间稳定时,将地面逆温的高度作为边界层高度。b风速极值法在风速廓线图上,可以自地面向上取风速达到地转风风速时的高度;或取达最大值的高度作为边界层高度。这里选用的是计算以风速达最大值时的高度作为边界层高度。1.4 模式方案设计 本实验方案的中心位于SACOL站(35.95N,104.13E),对A盆地地区的边界层结构进行模拟,如图1。图1.1 模拟区域地形高

9、度图 水平方向无嵌套,模拟区域为371*148,网格距为3公里。垂直方向采用地形追随坐标系,总共分46层,模式顶气压为10hPa,其中2000米以下分为19层。采用NCEP11再分析资料作为初始资料,每六小时输入一次,地形资料为NCEP提供的全球30s地形数据,每一小时输出一次模拟结果。WRF模式主要物理参数化过程包括WSM5微物理方案,RRTM长波辐射方案,Dudhia短波辐射方案,近地面层采用Monin-Obukhov方案,陆面过程采用Noah方案,边界层参数化方案主要选取YSU和MYJ方案两种,两种方案简介如下:在WRF模式的边界层参数化方案中,次网格尺度的湍流通量用诊断平均变量C通过垂

10、直扩散方程来参数化,即K理论闭合(一阶闭合理论): (3)其中,是变量C的湍流交换系数。aYSU方案YSU方案是在K理论的基础上建立起来的,因此被成为一阶闭合方案。在自由大气中,表示局地风切变和局地理查逊数的函数;在对流边界层中,是基于相似理论得到的。YSU方案除了考虑局地扩散外,还考虑了对流大尺度湍涡导致的非局地混合作用以及边界层顶的夹卷过程,即: (4)bMYJ方案MYJ方案需要TKE的诊断方程来使K理论方程闭合求解,因此该方案又被称为TKE闭合方案,也叫做湍流动能的1.5阶闭合方案。在这两种方案中,湍流交换系数通常被表述为 (5)其中:是混合长长度,是比例系数。该方案只考虑了局地混合作用

11、,不论是在对流边界层还是稳定边界层,垂直方向整层大气都采用上式计算湍流交换系数。【21】边界层结构特征的模拟分析边界层高度,边界层内位温廓线,比湿廓线,感热通量和潜热通量等的变化特征。接下来,选取2019年7月56日和2019年1月56日分别作为夏季和冬季的个例,对边界层的特征进行数值模拟。在模拟期间,天气晴朗,没有大型天气系统的存在与过境,不需要考虑云和降水的影响。此外,所选取的个例比较全面的包含了7:00和19:00的稳定度类型,具有很好的代表性。表1.1 模拟时段大气稳定度统计表时间稳定度类型时间稳定度类型2019.7.5-7:00B2019.1.5-7:00D2019.7.519:00

12、DE2019.1.519:00D2019.7.67:00D2019.1.67:00C2019.7.619:00D2019.1.619:00E1.1 边界层高度由表1.2可以看出,WRF模式可以较好的模拟出边界层的高度,模拟结果具有一定的可信度。在8次模拟实验中,YSU方案的模拟效果比较理想,在夏季7月6日7:00和19:00的模拟值与实测值的相对误差分别为0.9%和2.2%,在冬季1月5日7:00和19:00的模拟值与实测值的相对误差分别为12.1%和8.7%,可见,YSU方案对中性状态下的边界层高度的模拟效果非常好,这和YSU方案计算边界层高度的方法有关,YSU方案将第一个达到中性层结的高度

13、定为边界层高度,此外,YSU方案对夏季的模拟效果优于冬季。而MYJ方案的模拟结果则比较差,除7月5日7:00、 7月6日7:00和7月6日19:00三次模拟结果误差在20%之内以外,其它5次实验的模拟结果与实测值的偏差太大, 这可能是由于MYJ方案本身所考虑的物理过程十分复杂,受许多物理参数化方案影响,而本实验除边界层参数化方案单一变量外,其它方案均已设定好,由此导致许多假设条件不满足,从而影响结果。表1.2 两种边界层参数化方案模拟值与实测值对比表时间参数化方案模拟值实测值相对误差时间参数化方案模拟值实测值相对误差7.5-7:00YSUMYJOBS231.2244.6208.011.2%17

14、.6%1.5-7:00YSUMYJOBS80.021.991.012.1%75.9%7.5-19:00YSUMYJOBS554.5125.2481.015.3%74.0%1.5-19:00YSUMYJOBS94.022.0101.08.7%78.6%7.6-7:00YSUMYJOBS266.5246.3269.00.9%8.4%1.6-7:00YSUMYJOBS187.8110.1212.011.4%48.1%7.6-19:00YSUMYJOBS891.2756.5911.02.2%17.0%1.6-19:00YSUMYJOBS108.6438.6152.028.6%188.6%1.2 位温模

15、拟1.2.1 夏季如图1.1,a、b、c、d分别为2019年7月5日7:00、19:00和7月6日7:00、19:00的实测位温和模拟位温廓线图。在a图中,2019年7月5日7:00的边界层高度为208m,此时大气状态为中等不稳定状态,边界层内位温廓线属于一次线型。由图a可以看出,两种方案均模拟出了位温实测廓线的线性特征,在近边界层顶处,MYJ方案的模拟廓线出现了增长的趋势,与实测值偏离,相比较之下,YSU方案的模拟廓线与实测廓线更吻合,这是因为YSU方案中包含了边界层的夹卷作用,因此与实际情况更接近。在表1.3中的均方根误差值同样证明了这一点,YSU方案模拟值与实测值的均方根误差为0.168

16、6,而MYJ方案的为0.4682,可见YSU方案模拟效果更好一些。图1.1 夏季两种边界层参数化方案的位温模拟廓线与实测廓线对比图(a)2019年7月5日7:00;(b)2019年7月5日19:00;(c)2019年7月6日7:00;(d)2019年7月6日19:00;在b图中,2019年7月5日19:00的边界层高度为481m,此时的大气状况为微弱稳定型(DE),边界层内位温廓线型为逆温突变型。在边界层高度范围内,两种方案均可较好的模拟出位温的垂直分布。此时,在近地层100m一下,存在逆温,MYJ方案对逆温的模拟比YSU方案要好,主要是因为MYJ方案的局地闭合方法适合于研究近地面逆温的稳定状

17、态,但是在数值上,MYJ方案的模拟值与实测值相差较大,大约在1K左右。在100米以上至边界层顶范围内,YSU方案的模拟结果与实测值比较一致,有对流存在,上下位温比较均匀。表1.3中显示,YSU方案的模拟值与实测值的均方根误差为0.1844,而MYJ方案的误差为0.7692。相比较之下,YSU方案明显更优。在c图中,2019年7月6日7:00的边界层的高度为269m,此时大气的状况为中性(D),边界层内的位温廓线型属于逆温突变型。在地面至30米范围内,有逆温存在,YSU方案和MYJ方案均模拟出了逆温的存在。在30米以上,至边界层顶,存在较弱的对流,此时,太阳刚刚出来,太阳辐射比较弱,大气层结为中

18、性,并向不稳定发展,YSU方案在这一段模拟效果比较好,这主要是因为非局地闭合方案湍流交换较局地闭合方案能够较好的模拟对流边界层特点。在表1.3中,YSU方案的模拟值与实测值的均方根误差为0.3752,MYJ方案的模拟误差为0.4906,相比较之下,YSU方案的模拟结果更优。在d图中,2019年7月6日19:00的边界层高度为911m,此时大气状态为中性(D),大气处于对流不稳定向稳定状态的过渡时期,边界层内的位温廓线型为逆温突变型。自地面至大约离地100m高度的范围内,存在逆温,由图中,YSU方案的模拟廓线在边界层内为无逆温型,MYJ方案的模拟廓线在近地面有微弱的逆温,这可能是由于两种方案的响

19、应滞后造成的。表1.3中,YSU方案的模拟值与实测值的均方根误差为1.1011,MYJ方案的模拟误差为0.5080,因此,MYJ方案的模拟效果较好一些。表1.3 夏季两种边界层参数化方案位温模拟值与实测值的均方根误差YSUMYJ7.5-7:007.5-19:00 RMSE RMSE0.16860.18440.46820.76927.6-7:00 7.6-19:00 RMSE RMSE0.37521.10110.49060.50801.2.2 冬季如图1.2,a、b分别为2019年1月5日7:00、19:00和1月6日7:00、19:00的实测位温和模拟位温廓线图。在a图中,2019年1月5日7

20、:00的边界层高度为91m,此时大气状态为中性(D),边界层内的位温廓线为无逆温型。自地面至边界层顶,逆温已经被破坏,对流发展,两种方案均模拟出了无逆温型的特点,垂直方向上下位温混合均匀。但是在两种模拟结果与实测廓线均有一定的偏差,模拟效果一般。表1.4中,两种模拟结果与实测值的均方根误差分别为1.5617和1.6327,同样可以看出模拟效果不是很理想。在b图中,2019年1月5日19:00的边界层高度为103m,大气状态为中性(D),边界层内的位温廓线为逆温突变型。在近地面,存在很弱的逆温,这说明,地面的辐射冷却刚刚开始发展。YSU方案模拟出了逆温,而MYJ方案模拟出了一个非常强的逆温,与实

21、际情况不符。在逆温层顶至边界层顶这一段内,对流仍然存在。表1.4中,YSU方案和MYJ方案的模拟结果较差,分别为2.6402和4.8894。这个问题,有待进一步改进。图1.2 冬季两种边界层参数化方案的位温模拟廓线与实测廓线对比图(a)2019年1月5日7:00;(b)2019年1月5日19:00;(c)2019年1月6日7:00;(d)2019年1月6日19:00;在c图中,2019年1月6日7:00的边界层高度为212m,大气状态为弱不稳定(C),边界层内的位温廓线为逆温突变型。在地面至80米范围内,有逆温存在,而从逆温层顶至边界层顶,对流开始发展。YSU方案没有模拟出逆温的存在,在近地面

22、,呈现出混合层结构,位温垂直分布均匀,而MYJ方案模拟出较强的位温,且逆温层高度较高。总体看来,MYJ方案的模拟效果很好。在表1.4中,YSU方案的模拟误差为1.0385,MYJ方案的模拟误差为0.3490,MYJ方案的误差较小,主要是因为大气状态比较稳定,局地闭合方案的MYJ方案对位温的模拟效果好一些。在d图中,2019年1月6日19:00的边界层高度为152m,大气状态为弱稳定(E),边界层内的位温廓线为梯度渐变减弱型。在边界层高度范围内,位温梯度随高度减小,位温随高度的升高,增长越来越慢。两种方案在近地面模拟出了逆温的存在,但是逆温存在的高度仅为50m一下,结合表1.4,可以看出,两种方

23、案的模拟值与实测值的均方根误差均大于2,模拟效果也比较差。表1.4 冬季两种边界层参数化方案位温模拟值与实测值的均方根误差YSUMYJ 1.5-7:001.5-19:00 RMSE RMSE1.56172.64021.63274.88941.6-7:001.6-19:00 RMSE RMSE1.03852.33690.34902.54511.3 比湿模拟1.1.1 夏季如图1.3,a、b分别为2019年7月5日7:00、19:00和7月6日7:00、19:00的实测比湿和模拟比湿廓线图。图1.3 夏季两种边界层参数化方案的比湿模拟廓线与实测廓线对比图(a)2019年7月5日7:00;(b)20

24、19年7月5日19:00;(c)2019年7月6日7:00;(d)2019年7月6日19:00;在图a中,2019年7月5日7:00的边界层高度为208m,大气稳定度状态为中等不稳定(B),在边界层高度范围内,比湿廓线属于比湿突变型。从图中可以看出,两种方案的模拟值与实测值的模拟误差较大,表1.5中,两种方案模拟值与实测值的均方根误差均在2.6左右,模拟结果比较一般。在图b中,2019年7月5日19:00的边界层高度为481m,大气稳定度为微弱稳定(DE),在边界层高度范围内,比湿廓线型同样为比湿突变型。由图看出,两种方案模拟效果一般。实测比湿廓线呈现出先减小后增加的趋势,而YSU方案的模拟效

25、果体现为上下均匀,表现出了对流的特征。在从地面至200m高度范围内,比湿随高度递减,两种参数化方案中的MYJ方案模拟出了这种变化,但是在200m以上至边界层顶这一段,实际比湿增加,而MYJ方案的模拟值仍然随高度递减。在表1.5中,从数值上,两种方案的模拟值与实测值的均方根误差分别为1.2142和2.9999,模拟效果较差。在图c中,2019年7月6日7:00的边界层高度为269m,大气为中性(D),在边界层范围内,比湿廓线为比湿突变型。从实测廓线分析,边界层内比湿随高度递减,在40m出现突变,比湿递减率减小,从总体看,比湿随高度呈现出递减的特征,YSU方案和MYJ方案均模拟出了这种趋势。在表1

26、.5中,两种方案的均方根误差分别为0.4190和0.91994,从数值上看,YSU方案的模拟效果更优。在图d中,2019年7月6日19:00的边界层高度为911m,大气状态为中性(D),在边界层范围内,比湿递减率随高度减小,廓线型呈现为梯度渐变减弱型。两种参数化方案在边界层高度内模拟的效果表现为比湿随高度不变,与实测廓线偏差较大。在表1.5中,两种方案模拟的均方根误差在1.3左右,误差较大。表1.5 夏季两种边界层参数化方案比湿模拟值与实测值的均方根误差YSUMYJ7.5-7:007.5 19:00 RMSE RMSE2.63761.21422.67002.99997.6 7:007.6-19

27、:00 RMSE RMSE0.41901.34060.91941.29411.1.2 冬季如图1.4,a、b分别为2019年1月5日7:00、19:00和1月6日7:00、19:00的实测比湿和模拟比湿廓线图。在图a中,2019年1月5日7:00的边界层高度为91m,大气状态为中性(D),此时,边界层内实测比湿廓线型为一次线型。由图可以看出,两种参数化方案均较好的模拟出了比湿的变化趋势,在表1.6中,两种方案的模拟值与实测值的均方根误差分别为0.0283和0.0586,模拟效果非常好。图1.4 冬季两种边界层参数化方案的比湿模拟廓线与实测廓线对比图(a)2019年1月5日早晨7:00;(b)2

28、019年1月5日晚上19:00;(c)2019年1月6日早晨7:00;(d)2019年1月6日晚上19:00;在图b中,2019年1月5日19:00的边界层高度为103m,大气状态为中性(D),此时,边界层内的实测比湿廓线型为逆比湿突变型。在近地层,有逆比湿存在,在大约50m高度以上至边界层顶,比湿呈现出先减小的趋势,两种参数化方案均未模拟出逆比湿的存在,这可能是由于模式对于地面辐射冷却的响应不够敏捷,进而导致模拟出现了滞后的现象。在表1.6中,两种方案的模拟值与实测值的均方根误差分别为0.3948和0.3934,从数值上看模拟效果较好。在图c中,2019年1月6日7:00的边界层高度为212

29、m,大气状态为弱不稳定(C),此时,边界层内的实测比湿廓线型为逆比湿突变型。从实测数据廓线看,近地层存在逆比湿,且在30m高度往上,比湿随高度开始降低。YSU方案和MYJ方案同样未模拟出了这种趋势。在表1.6中,两种方案的均方根误差分别为0.2305和0.1354,从数值上看,MYJ方案的模拟效果更优。在图d中,2019年1月6日19:00的边界层高度为152m,大气状态为弱稳定(E),此时,边界层内的实测比湿廓线型为逆比湿突变型。在近地面50m高度以下,存在逆比湿,且在此高度往上,比湿随高度减小。从图中可以看出,两种参数化方案均模拟出了逆比湿的存在,其中,MYJ方案模拟的逆比湿高度较高,且逆

30、比湿的递增率较实测值大。而YSU方案的模拟效果与实测廓线更为接近。从表1.6中,可以看出两种方案的模拟误差均在0.3附近,并且,YSU方案的模拟效果相对更优一些。表1.6 冬季两种边界层参数化方案比湿模拟值与实测值的均方根误差YSUMYJ1.5-7:001.5-19:00RMSERMSE0.02830.39480.05860.39341.6-7:001.6-19:00RMSERMSE0.23050.31780.13540.34961.4 热通量模拟1.4.1 感热通量如图1.5,是夏季和冬季两种边界层方案YSU方案和MYJ方案模拟的SACOL近地面感热通量与实测资料的对比情况。图1.5 7月5

31、6日(a)和1月56日(b)两种边界层参数化方案的感热通量模拟值与实测值对比图感热是湍流运动从地面向大气传输热量通量的形式,其受下垫面影响很大。由图1.5(a)、(b)可以看出,WRF模式模拟的感热通量的变化趋势与实测值基本吻合,对低值的模拟效果非常好,但是在两个峰值处的模拟有一定的误差,其模拟结果与近地层的温度的模拟结果相似,感热的计算主要是依赖近地面的温度和地表的植被覆盖状况,故对温度的模拟结果与对感热通量的模拟结果相关。夏季日出之后,随着太阳高度角的增大,感热通量逐渐增大,在7月5日6:00之后出现了正值,此时,能量是由地面向大气传输,在中午13:00,感热通量达到最大值,为328.6W/,此后,在21:

copyright@ 2008-2023 冰点文库 网站版权所有

经营许可证编号:鄂ICP备19020893号-2