天气学诊断实习六 计算水汽流函数和水汽势函数

天气学诊断实习六 计算水汽流函数和水汽势函数 一、实习目的:

熟悉水汽流函数和水汽势函数程序的实际编程计算。

二、实习内容:

编制计算水汽流函数和水汽势函数的程序，并且绘制出两个时次25日20时，26日20时的水汽流函数和水汽势函数分布（850hPa,500hPa)

三、算法原理

势函数初值： 0
流函数初值：

利布曼法：

四、代码实现

# -*- 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 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 hd(x):
    a=math.pi/180*x
    return a

def sd(u,v,leftlon, rightlon, lowerlat, upperlat,f):
    a=6371000
    b=hd(f)
    n1=len(u[:,0])-1
    n2=len(u[0,:])-1
    N=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 wd(u,v,leftlon, rightlon, lowerlat, upperlat,f):
    a=6371000
    b=hd(f)
    n1=len(u[:,0])-1
    n2=len(u[0,:])-1
    N=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)*((v[1,j]-v[0,j])/(math.cos(hd(lat[0]))*b)-(u[0,j+1]-u[0,j-1])/b+2*u[0,j]*math.tan(hd(lat[0])))
            data[n1,j]=(1/2/a)*((v[n1-1,j]-v[n1,j])/(math.cos(hd(lat[n1]))*b)-(u[n1,j+1]-u[n1,j-1])/b+2*u[n1,j]*math.tan(hd(lat[n1])))
            data[i,0]=(1/2/a)*((v[i+1,0]-v[i-1,0])/(math.cos(hd(lat[i]))*b)-(u[i,1]-u[i,0])/b+2*u[i,0]*math.tan(hd(lat[i])))
            data[i,n2]=(1/2/a)*((v[i+1,n2]-v[i-1,n2])/(math.cos(hd(lat[i]))*b)-(u[i,n2]-u[i,n2-1])/b+2*u[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)*((v[i+1,j]-v[i-1,j])/(math.cos(hd(lat[i]))*b)-(u[i,j+1]-u[i,j-1])/b+2*u[i,j]*math.tan(hd(lat[i])))
    ##四角差分
    data[0,0]=(1/2/a)*((v[1,0]-v[0,0])/(math.cos(hd(lat[0]))*b)-(u[0,1]-u[0,0])/b+2*u[0,0]*math.tan(hd(lat[0])))
    data[0,n2]=(1/2/a)*((v[1,n2]-v[0,n2])/(math.cos(hd(lat[0]))*b)-(u[0,n2-1]-u[0,n2])/b+2*u[0,n2]*math.tan(hd(lat[0])))
    data[n1,0]=(1/2/a)*((v[n1-1,0]-v[n1,0])/(math.cos(hd(lat[n1]))*b)-(u[n1,0]-u[n1,1])/b+2*u[n1,0]*math.tan(hd(lat[n1])))
    data[n1,n2]=(1/2/a)*((v[n1-1,n2]-v[n1,n2])/(math.cos(hd(lat[n1]))*b)-(u[n1,n2-1]-u[n1,n2])/b+2*u[n1,n1]*math.tan(hd(lat[n1])))
    return data

def shi(D,leftlon, rightlon, lowerlat, upperlat,f):
    kapa=np.zeros((29,53))
    R=np.zeros((29,53))
    a=6371000
    b=hd(f)
    N=int((upperlat-lowerlat)/f+1)
    lat=np.linspace(lowerlat, upperlat,N)
    lat= lat[::-1]
    for k in range(0,500):                  ##迭代次数
        t=kapa[1,1]
        #中间差分迭代部分
        for i in range(1,28):
            for j in range(1,52):
                R[i,j]=(kapa[i+1,j]+kapa[i-1,j])/((a*math.cos(hd(lat[i]))*b)**2)+(kapa[i,j+1]+kapa[i,j-1])/((a*b)**2)+D[i,j]-(2/((a*math.cos(hd(lat[i]))*b)**2)+2/(a*b)**2)*kapa[i,j]
                kapa[i,j]=kapa[i,j]+R[i,j]/2/(1/((a*math.cos(hd(lat[i]))*b)**2)+1/(a*b)**2)             
     
        if (k>1) & (R[1,1]1) & (R[1,1]1) & (R[1,1]1) & (R[1,1] 
五、实习结果 

 
欢迎分享，转载请注明来源：内存溢出
原文地址:https://54852.com/zaji/5651180.html