计算假相当位温

wsx

本人计算假想当位温的程序，公式已经核对过了:begin;=================================================================================================
f = addfile("data.nc", "r")
time = f->time
ym = cd_calendar(time,-3)

hour = ym%100
day  = ym/100%100
hour = hour+8
day1  = where(hour.ge.24,day+1,day)    ;
hour1 = where(hour.ge.24,hour-24,hour) ;
idx = ind(day1.ge.17.and.day1.le.24)   ;转为北京时，并选取时间

;==========================假相当位温计算过程===========================================================
z = short2flt(f->z(idx,:,:,:))
temp = dimsizes(z)
se2 = new((/temp(0),temp(1),temp(2),temp(3)/),float)
pres = int2flt(f->level );;;;;“重点注意”
do i = 0,dimsizes(pres)-1
prs = pres(i)
tk  = short2flt(f->t(idx,{prs},:,:));温度 ;K
rh  = short2flt(f->r(idx,{prs},:,:));相对湿度;%
tc = tk - 273.16  ;转换为摄氏温度
es =6.112*exp(17.67*tc/(tc+243.5)) ;修正后的Tetens经验公式,单位hPa;  饱和水汽压
;es = where( tc.lt.0, 6.112*exp(17.67* tc /( tc +243.5)),\
;                           6.11*10^(7.5* tc /( tc +237.3)) )
qs = 0.622*es/(prs - 0.378*es) ;计算饱和比湿，单位g/g   
;e = es * rh/100  ;计算水汽压，根据rh的单位原因除以100   相对湿度 = 水汽压/饱和水汽压
q  = short2flt(f->q(idx,{prs},:,:));比湿
e     =  pres(i)*q/(0.622+q)
tlcl = 55.0+2840.0/(3.5*log(tk)-log(e)-4.805) ;凝结高度的绝对温度
theta = tk*((1000/prs)^(0.2854*(1.0-0.28*q)))
se = theta*exp(((3376./tlcl)-2.54)*q*(1.0+0.81*q))-273.16;假相当位温
;================================================================================================
se2(:,i,:,:) = se
end do
copy_VarMeta(z,se2)

;================================================================================================
lat_station = 37.77
lon_station = 106.16
se3 = se2(:,{1000:200},{lat_station},{lon_station});time , level


;================================================================================================
wks = gsn_open_wks("pdf","se_day20")
rh_res              = True
  rh_res@cnSmoothingOn  = True
  rh_res@cnSmoothingDistanceF = 0.005
  rh_res@cnSmoothingTensionF = 0.01
  rh_res@trYReverse   = True     ; Reverse the Y values.;第一层是200hPa,所以要用这个语句倒转Y轴
  ;rh_res@trXReverse   = True
  rh_res@gsnDraw      = False    ; Don't draw individual plot  不画单独的图层
  rh_res@gsnFrame     = False    ; Don't advance frame.
  rh_res@vpXF         = 0.1    ; x location   X轴左边界的位置
  ;rh_res@vpYF         = 0.90     ; y location   Y轴上边界的位置
  rh_res@vpWidthF     = 0.83      ; width       图片的宽度
  rh_res@vpHeightF    = 0.40     ; height
  rh_res@tmXBLabelFontHeightF = 0.015   ;设置xy轴字体大小
  rh_res@tmYLLabelFontHeightF = 0.015
  rh_res@cnFillOn          = True               ; turns on color fill
  rh_res@cnLinesOn         = False
  rh_res@gsnLeftString =" "
  rh_res@gsnLeftStringFontHeightF =0.015
  rh_res@gsnRightString =" "
  rh_res@gsnRightStringFontHeightF =0.015
  rh_res@cnLevelSelectionMode       = "ExplicitLevels"   ;AutomaticLevels
  rh_res@cnLevels=(/-5,0,5,10,15,20,25,30,35,40,45,50,55,60,65,70/)
  rh_res@cnLineLabelsOn    = False               ; no contour labels
  rh_res@tiMainString      = " " ; title
  rh_res@tiMainFontHeightF    = 0.02
  rh_res@lbLabelFontHeightF      =  0.01                   ;-- increase label font size

  rh_res@tmYLMode       = "Explicit"
  rh_res@tmYLLabels     = (/"300", "400", "500", "600","700","800", \
                            "850","900","925","1000"/)

  time2 = se3&time
;---------------------------------------------------------------------------------------------------------------
  rh_res@tmXBMode = "Explicit"
  time_x = new(31,float)
  do i = 0,30
  time_x(i) = time2(i*6)
  end do
  rh_res@tmXBValues = (/time_x/)
  rh_res@tmXBLabels     = (/"8", "14", "20", "2","8", "14", "20", "2","8", "14", "20", "2","8", "14", "20", "2",\
  "8", "14", "20", "2","8", "14", "20", "2","8", "14", "20", "2","8", "14","20"/)

  rh_res@tmXTOn   = False
  rh_res@tmYROn   = False



  rh_res@gsnXRefLine           = (/time2(16),time2(40),time2(64),time2(80),time2(112),time2(136),time2(160)/)            ; create a reference line   
  rh_res@gsnYRefLineThicknesses = 5
  rhfill    = gsn_csm_contour(wks,se3(level|:,time|:),rh_res)
  draw(rhfill)               
  frame(wks)
end

jinwithsky

谢谢分享~刚好需要