- 积分
- 299
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-4-6
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
; These files are loaded by default in NCL V6.2.0 and newer
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
FILE = systemfunc("ls"+" /cygdrive/L/DJBY/20y_winter/2018_1/*.grib2")
numFILE = dimsizes(FILE)
mpv2z= new((/181,360/),float,0.)
mpv1z= new((/181,360/),float,0.)
do j = 6,10
file_list = addfile(FILE(j), "r")
time=str_get_cols(FILE(j),39,49)
lon = file_list->lon_0
lat = file_list->lat_0
lon@units="degrees_east"
lat@units="degrees north"
dlon = (lon(2)-lon(1))*0.0174533
dlat = (lat(2)-lat(1)) * 0.0174533
;print(time)
nlev1=750
tmpx = file_list->TMP_P0_L100_GLL0({nlev1*100},:,:)
rhx = file_list->RH_P0_L100_GLL0({nlev1*100},:,:)
prsx=nlev1
esx=(6.112*exp(17.67*(tmpx-273.15)/(tmpx-29.65)))
qx=rhx*(0.62197*esx/(prsx-esx))/100.
ex=prsx*qx/(0.62197+qx)+1e-10
tlclx=55.0+2840.0/(3.5*log(tmpx)-log(ex)-4.805)
thetax=tmpx*((1000/prsx)^(0.2854*(1.0-0.28*qx)))
eqtx=thetax*exp(((3376./tlclx)-2.54)*qx*(1.0+0.81*qx))
vx= file_list->VGRD_P0_L100_GLL0({nlev1*100},:,:)
ux= file_list->UGRD_P0_L100_GLL0({nlev1*100},:,:)
nlev2 = 650
tmpy = file_list->TMP_P0_L100_GLL0({nlev2*100},:,:)
rhy = file_list->RH_P0_L100_GLL0({nlev2*100},:,:)
prsy=nlev2
esy=(6.112*exp(17.67*(tmpy-273.15)/(tmpy-29.65)))
qy=rhy*(0.62197*esy/(prsy-esy))/100
ey=prsy*qy/(0.62197+qy)+1e-10
tlcly=55.0+2840.0/(3.5*log(tmpy)-log(ey)-4.805)
thetay=tmpy*((1000/prsy)^(0.2854*(1.0-0.28*qy)))
eqty=thetay*exp(((3376./tlcly)-2.54)*qy*(1.0+0.81*qy))
vy = file_list->VGRD_P0_L100_GLL0({nlev2*100},:,:)
uy = file_list->UGRD_P0_L100_GLL0({nlev2*100},:,:)
nlev3 = 700
tmpz = file_list->TMP_P0_L100_GLL0({nlev3*100},:,:)
rhz = file_list->RH_P0_L100_GLL0({nlev3*100},:,:)
absvprs = file_list->ABSV_P0_L100_GLL0({nlev3*100},:,:)
prsz = nlev3
esz = (6.112*exp(17.67*(tmpz-273.15)/(tmpz-29.65)))
qz = rhz*(0.62197*esz/(prsz-esz))/100.
ez = prsz*qz/(0.62197+qz)+1e-10
tlclz = 55.0+2840.0/(3.5*log(tmpz)-log(ez)-4.805)
thetaz = tmpz*((1000/prsz)^(0.2854*(1.0-0.28*qz)))
eqtz = thetaz*exp(((3376./tlclz)-2.54)*qz*(1.0+0.81*qz))
rhz@_FillValue = -999.0
numlat = dimsizes(lat)
numlon = dimsizes(lon)
dtdX=new((/numlat, numlon/), typeof(rhz), rhz@_FillValue)
dtdY=new((/numlat, numlon/), typeof(rhz), rhz@_FillValue)
do n1=0,numlat-1
dX = 6378388.*cos(0.0174533*lat(n1))*dlon ; constant at this latitude
dtdX(n1:n1,:) = center_finite_diff_n(eqtz(n1:n1,:), dX ,False,0,1)
end do
do n2=0,numlon-1
dY = 6378388. * dlat
dtdY(:,n2:n2) = center_finite_diff_n(eqtz(:,n2:n2), dY, False, 0, 0)
end do
eqtxy=eqtx-eqty
vxy=vx-vy
uxy=ux-uy
dp=10000
mp1=-9.8*absvprs*eqtxy/dp
mp2=9.8*((vxy/dp)*dtdX-(uxy/dp)*dtdY)
mpv=mp1+mp2
mpv2=mp2*1E6
mpv1=mp1*1E6
mpv2z=(mpv2z+mpv2)/5.
mpv1z=(mpv1z+mpv1)/5.
end do
copy_VarCoords(rhz,mpv1z)
copy_VarCoords(rhz,mpv2z)
wks = gsn_open_wks("png","L:/DJBY/ca_saf/201801/MPV12/mpv12z_"+time)
hgtres = True
hgtres@cnSmoothingOn=True
hgtres@cnSmoothingDistanceF=0.0001
hgtres@cnSmoothingTensionF=-0.2
hgtres@mpMaxLatF = 30 ; specify the plot domain
hgtres@mpMinLatF = 20. ;
hgtres@mpMinLonF = 96. ;
hgtres@mpMaxLonF = 108. ;
hgtres@pmTickMarkDisplayMode = "Never"
hgtres@gsnMajorLonSpacing =10.
hgtres@gsnMinorLonSpacing =2.
hgtres@gsnMajorLatSpacing =10
hgtres@gsnMinorLatSpacing =2.
hgtres@mpOutlineOn = True ; turn the map outline on
hgtres@gsnLeftString = ""
hgtres@gsnRightString = ""
hgtres@gsnDraw = False ; don't draw
hgtres@gsnFrame = False ; don't advance frame
hgtres@mpDataSetName = "Earth..4"
hgtres@mpDataBaseVersion = "MediumRes"
hgtres@mpOutlineOn = True
hgtres@mpOutlineSpecifiers = (/"China:states","Taiwan"/)
hgtres@mpGeophysicalLineThicknessF = 2
hgtres@mpNationalLineThicknessF = 2
hgtres@cnLinesOn = False ; turn of contour lines
hgtres@cnFillPalette = "BlueWhiteOrangeRed"
hgtres@cnFillOn = True
hgtres@cnInfoLabelOn = False ; disappear lable information
hgtres@lbBoxMinorExtentF =0.2
hgtres@gsnLeftStringFontHeightF =0.015
hgtres@lbTitleFontHeightF =0.012
hgtres@lbLabelFontHeightF =0.013
hgtres@lbLabelStride =2
hgtres@cnLevelSelectionMode = "ManualLevels" ; manual contour levels
hgtres@cnMinLevelValF = -1.2 ; minimum level
hgtres@cnMaxLevelValF = 1.2 ; maximum level
hgtres@cnLevelSpacingF = 0.2 ; contour spacing
hgtres@tmXBTickStartF=96
hgtres@tmXBTickEndF=108
hgtres@tmXBTickSpacingF=2
hgtres@tmYLTickStartF=20
hgtres@tmYLTickEndF=30
hgtres@tmYLTickSpacingF=2
resc = True
resc@cnSmoothingOn=True
resc@cnSmoothingDistanceF=0.0001
resc@cnSmoothingTensionF=-0.2
resc@cnLinesOn = True
;resc@cnLevelSelectionMode ="ManualLevels"
resc@cnMaxLevelValF =0.1
resc@cnMinLevelValF =-0.1
resc@cnLevelSpacingF = 0.05
resc@cnLineThicknessF =2.
resc@cnLineColor ="black"
resc@gsnContourNegLineDashPattern = 11
resc@gsnContourPosLineDashPattern =0
plot_mpv1 = gsn_csm_contour_map(wks,mpv1z,hgtres)
plot_mpv2 = gsn_csm_contour(wks,mpv2z,resc) ; create the U-wind plot
overlay(plot_mpv1,plot_mpv2) ; overlay the U-wind plot on the temperature plot
draw(plot_mpv1) ; draw the temperature plot (with the U-wind plot overlaid)
frame(wks) ; advance the frame
end
|
-
|