天气学诊断实习四计算垂直速度

一、实习目的:

熟悉垂直差分在气象中的应用，掌握垂直速度的实际编程计算。

二、实习内容:

编制计算垂直速度程序，并绘制500hPa垂直速度。用第二种修正方案，其中大气层顶的垂直速度可以直接采用0，也可以用绝热法。并且绘制出两个时次25日20时，26日20时的修正后的垂直速度分布（850hPa,500hPa)

三、算法原理:

垂直速度计算：

修正方案二：

四、代码实现

# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd
import math
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
from pylab import *                                 #支持中文
mpl.rcParams['font.sans-serif'] = ['SimHei']##角度转弧度
def hd(x):a=math.pi/180*xreturn a
##micaps读数据函数
def micaps(a):data=np.zeros((29,53))for i in range(0,29):data[i,0:10]=a[i*6,0:10]data[i,10:20]=a[i*6+1,0:10]data[i,20:30]=a[i*6+2,0:10]data[i,30:40]=a[i*6+3,0:10]data[i,40:50]=a[i*6+4,0:10]data[i,50:53]=a[i*6+5,0:3]return data
##散度计算
def sd(u,v,leftlon, rightlon, lowerlat, upperlat,f):a=6371000b=hd(f)n1=len(u[:,0])-1n2=len(u[0,:])-1N=int((upperlat-lowerlat)/f+1)lat=np.linspace(lowerlat, upperlat,N)lat= lat[::-1]data=np.zeros((29,53))##四边差分for i in range(1,n1):for j in range (1,n2):data[0,j]=(1/2/a)*((u[1,j+1]-u[0,j-1])/(math.cos(hd(lat[0]))*b)+(v[1,j]-v[0,j])/b-2*v[0,j]*math.tan(hd(lat[0])))data[n1,j]=(1/2/a)*((u[n1,j+1]-u[n1,j-1])/(math.cos(hd(lat[n1]))*b)+(v[n1-1,j]-v[n1,j])/b-2*v[n1,j]*math.tan(hd(lat[n1])))data[i,0]=(1/2/a)*((u[i,1]-u[i,0])/(math.cos(hd(lat[i]))*b)+(v[i+1,1]-v[i-1,0])/b-2*v[i,0]*math.tan(hd(lat[i])))data[i,n2]=(1/2/a)*((u[i,n2-1]-u[i,n2])/(math.cos(hd(lat[i]))*b)+(v[i+1,n2]-v[i-1,n2])/b-2*v[i,n2]*math.tan(hd(lat[i])))##中间部分差分for i in range(1,n1):for j in range (1,n2):data[i,j]=(1/2/a)*((u[i,j+1]-u[i,j-1])/(math.cos(hd(lat[i]))*b)+(v[i+1,j]-v[i-1,j])/b-2*v[i,j]*math.tan(hd(lat[i])))##四角差分data[0,0]=(1/2/a)*((u[0,1]-u[0,0])/(math.cos(hd(lat[0]))*b)+(v[1,0]-v[0,0])/b-2*v[0,0]*math.tan(hd(lat[0])))data[0,n2]=(1/2/a)*((u[0,n2-1]-u[0,n2])/(math.cos(hd(lat[0]))*b)+(v[0,n2-1]-v[0,n2])/b-2*v[0,n2]*math.tan(hd(lat[0])))data[n1,0]=(1/2/a)*((u[n1,1]-u[n1,0])/(math.cos(hd(lat[n1]))*b)+(v[n1-1,0]-v[n1,0])/b-2*v[n1,0]*math.tan(hd(lat[n1])))data[n1,n2]=(1/2/a)*((u[n1-1,n2]-u[n1,n2])/(math.cos(hd(lat[n1]))*b)+(v[n1-1,n2]-v[n1,n2])/b-2*v[n1,n1]*math.tan(hd(lat[n1])))return data
##画图函数
def draw(data,leftlon, rightlon, lowerlat, upperlat,a,name):lon=np.arange(leftlon, rightlon+0.01,a)lat=np.arange(lowerlat, upperlat+0.01,a)data=data[::-1, :]##坐标反转#建立画布proj = ccrs.PlateCarree() # 设置投影fig, f2_ax1 = plt.subplots(figsize=(15,15), subplot_kw=dict(projection=proj))leftlon, rightlon, lowerlat, upperlat = (leftlon, rightlon, lowerlat, upperlat)#绘制data1=f2_ax1.contourf(lon,lat,data)# data2=f2_ax1.contour(lon,lat,data,colors='k', linewidths=1, linestyles='solid',levels=np.linspace(-32,60,23))# plt.clabel(data2,fontsize=10,colors='r',fmt='%.2f')#在画布的绝对坐标建立子图f2_ax1.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())#海岸线，50m精度f2_ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))#湖泊数据f2_ax1.add_feature(cfeature.LAKES, alpha=0.5)#以下6条语句是定义地理坐标标签格式f2_ax1.set_xticks(np.arange(leftlon,rightlon+10,10), crs=ccrs.PlateCarree())f2_ax1.set_yticks(np.arange(lowerlat,upperlat+10,10), crs=ccrs.PlateCarree())lon_formatter = cticker.LongitudeFormatter()lat_formatter = cticker.LatitudeFormatter()f2_ax1.xaxis.set_major_formatter(lon_formatter)f2_ax1.yaxis.set_major_formatter(lat_formatter)f2_ax1.set_title(name,loc='center',fontsize =20)# shrink 控制 colorbar 长度，pad 控制colorbar和图的距离plt.rcParams['axes.unicode_minus'] = False##负号显示问题plt.colorbar(data1, shrink=0.3, pad=0.02)#orientation='horizontal'位置参数namespace = globals()
a=[1000,925,850,700,500,400,300,250,200,150,100]
#读数据
for x in a:namespace['uv_%d' % x] = pd.read_table(r'C:\Users\马冠龙\Desktop\micaps 11层\micaps 11层\uv\%d\13052620.000'%(x),header=None,skiprows=3,sep='\s+')namespace['uv_%d' % x] = np.array(namespace['uv_%d' % x])namespace['u_%d' % x]=micaps(namespace['uv_%d' % x][0:174,:])namespace['v_%d' % x]=micaps(namespace['uv_%d' % x][174:,:])#散度计算
for x in a:namespace['D%d' % x]=sd(namespace['u_%d' % x],namespace['v_%d' % x],30,160,10,80,2.5)#垂直速度计算
w_1000=np.zeros((29,53))
w_925=w_1000+0.5*(D1000+D925)*75
w_850=w_925+0.5*(D925+D850)*75
w_700=w_850+0.5*(D850+D700)*150
w_500=w_700+0.5*(D700+D500)*200
w_400=w_500+0.5*(D500+D400)*100
w_300=w_400+0.5*(D400+D300)*100
w_250=w_300+0.5*(D300+D250)*50
w_200=w_250+0.5*(D250+D200)*50
w_150=w_200+0.5*(D200+D150)*50
w_100=w_150+0.5*(D150+D100)*50
##垂直速度修正
M=0.5*11*10
k=2
w_850=w_850-k*(k+1)/2/M*(w_100-0)
k=4
w_500=w_500-k*(k+1)/2/M*(w_100-0)
draw(w_500,30,160,10,80,2.5,'500hPa垂直速度图  2013年5月26日20时')
draw(w_850,30,160,10,80,2.5,'850hPa垂直速度图  2013年5月26日20时')

五、实习结果

天气学诊断实习四计算垂直速度相关推荐

天气学诊断实习五计算水汽平流、水汽通量、水汽通量散度
一.实习目的: 熟悉水汽平流.水汽通量.水汽通量散度的实际编程计算. 二.实习内容: 编制计算水汽平流.水汽通量.水汽通量散度的程序,并且绘制出两个时次25日20时,26日20时的水汽平流.水汽通量. ...
天气学诊断实习六计算水汽流函数和水汽势函数
一.实习目的: 熟悉水汽流函数和水汽势函数程序的实际编程计算. 二.实习内容: 编制计算水汽流函数和水汽势函数的程序,并且绘制出两个时次25日20时,26日20时的水汽流函数和水汽势函数分布(850h ...
天气学诊断实习三计算涡度、散度、涡度平流和温度平流
一.实习目的: 熟悉使用差分方法解决方程中的导数.偏导数,计算涡度.散度.涡度平流和温度平流的程序二.实习内容: 编制计算涡度.散度.涡度平流和温度平流的程序,并且绘制出两个时次25日20时,26日 ...
天气学诊断实习一露点和凝结抬升高度的计算
一.实习内容已知南京站某一时次的相对湿度(92.8%).气温(14.1ºC).气压(1017.4, hPa)及站点高度(35.2m),请利用迭代法求出: (1) 该站的露点温度 (2) 凝结高度 ...
天气学诊断分析I 实习报告(实习四)
实习4 利用风场和高度场计算运动特征参量实习目的熟悉地转风涡度在气象中的应用,掌握涡度和散度的实际编程计算. 实习内容已知2020年的6h再分析资料,要素场有温度.高度.相对湿度和水平风场.求 ...
天气学诊断分析I 实习报告（二）
实习目的熟悉迭代方法在气象中的应用,掌握稳定度指数的实际编程计算. 实习内容 1.已知2020年的6h再分析资料,要素场有温度.高度.相对湿度和水平风场.用迭代法求某一区域(10ºN-60ºN,60 ...
天气学诊断分析I 实习报告（三）
实习目的熟悉迭代方法在气象中的应用,掌握稳定度指数的实际编程计算. 实习内容 1.已知2020年的6h再分析资料,要素场有温度.高度.相对湿度和水平风场.用迭代法求某一区域(10ºN-60ºN,60 ...
天气学原理和方法第四版pdf_天气学原理和方法汇总很好很全面.pdf
天气学原理和方法汇总很好很全面天气学原理和方法第一章大气运动的基本特征 1.大气运动受什么定律支配? 质量守衡.动量守衡和能量守衡定律 2.影响大气运动的真实力有哪几种? 气压梯度力.地心引力 ...
python计算矩形周长_一边学编程，一边学语数外，用python编程学三年级周长计算...
原标题:一边学编程,一边学语数外,用python编程学三年级周长计算编程并不神秘编程只是解决问题的一共方法 python是一门编程语言 python是一种解决问题的编程工具在小学阶段,学习编程的 ...

天气学诊断实习四计算垂直速度

一、实习目的:

二、实习内容:

三、算法原理:

四、代码实现

五、实习结果

天气学诊断实习四计算垂直速度相关推荐

最新文章

热门文章

天气学诊断实习四 计算垂直速度

一、实习目的:

二、实习内容:

三、算法原理:

四、代码实现

五、实习结果

天气学诊断实习四 计算垂直速度相关推荐

最新文章

热门文章

天气学诊断实习四计算垂直速度

天气学诊断实习四计算垂直速度相关推荐