ncl计算假相当位温和水汽通量

juliekxkl

看到ncl画假相当位温和水汽通量的帖子，和用grads画图脚本对比，发现有些不同。a1和b1的选取，where在此处为什么以263为界限选择不同公式，这个公式对吗？e1=P*q1/100./(0.62197+q1/100.0)后面需要叫上加上1e-10嘛？


load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin;
   a = addfile("./wrfout_d02_2006-06-05_00:00:00.nc","r"))  
;  type = "x11"
;  type = "pdf"
  type = "ps"
   wks = gsn_open_wks(type,"plt_PressureLevel3")
  gsn_define_colormap(wks,"wh-bl-gr-ye-re")
  res = True
   res@MainTitle                   = "REAL-TIME WRF"
   res@Footer = False

   pltres = True
   mpres = True
   mpres@mpGeophysicalLineColor      = "Black"
   mpres@mpNationalLineColor         = "Black"
   mpres@mpUSStateLineColor          = "Black"
   mpres@mpGridLineColor             = "Black"
   mpres@mpLimbLineColor             = "Black"
   mpres@mpPerimLineColor            = "Black"
   mpres@mpGeophysicalLineThicknessF = 2.0
   mpres@mpGridLineThicknessF        = 2.0
   mpres@mpLimbLineThicknessF        = 2.0
   mpres@mpNationalLineThicknessF    = 2.0
   mpres@mpUSStateLineThicknessF     = 2.0

   mpres@mpDataBaseVersion="MediumRes"
   mpres@mpDataSetName="Earth..4"
   mpres@mpOutlineSpecifiers=(/"China","Fujian","Guandong","Jiangxi","Zhejiang"/)
   times  = wrf_user_list_times(a)  ; get times in the filewww.mnmuc.org)
   ntimes = dimsizes(times)         ; number of times in the file

    tc = wrf_user_getvar(a,"tc",-1)        ; T in C
    u  = wrf_user_getvar(a,"ua",-1)        ; u averaged to mass points
    v  = wrf_user_getvar(a,"va",-1)        ; v averaged to mass points
    p  = wrf_user_getvar(a, "pressure",-1) ;
     rh = wrf_user_getvar(a,"rh",-1)        ; relative humidity;

   pressure_levels = (/ 1000.,950.,900.,850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,\
               350.,300.,250.,200.,150.,100./)   ; pressure levels to plot
   nlevels         = dimsizes(pressure_levels)     ; number of pressure levels
   tc_plane = wrf_user_intrp3d(tc,p,"h",pressure_levels,0.,False)
   z_plane  = wrf_user_intrp3d( z,p,"h",pressure_levels,0.,False)
   p  rh_plane = wrf_user_intrp3d(rh,p,"h",pressure_levels,0.,False)
   u_plane  = wrf_user_intrp3d( u,p,"h",pressure_levels,0.,False)
   p_plane  = wrf_user_intrp3d( p,p,"h",pressure_levels,0.,False)

   spd     = (u_plane*u_plane + v_plane*v_plane)^(0.5) ;
     spd@description = "Wind Speed"
     spd@units = "m/s"

;*计算假相当位温
     tk=tc_plane+273.166 k; W8 ]: G. [% v
     P = conform(tc_plane,pressure_levels ,1)
    print("p")
     a1 = where(tk .gt. 263.0, 0.622*6.11*exp(17.26*(tk-273.16)/(tk-35.86)), \
                     0.622*6.11*exp(21.87*(tk-273.16)/(tk-7.66)))   ;
     b1= where(tk .gt. 263.0, P-0.278*exp(17.26*(tk-273.16)/(tk-35.86)),\
                 P-0.278*exp(21.87*(tk-273.16)/(tk-7.66)))   ;
     qs1=a1/b1
q1=qs1*rh_plane
     e1=P*q1/100./(0.62197+q1/100.0)

     tk1=55.0+2840.0/(3.5*log(tk)-log(e1)-4.805)

     pot1=tk*(1000/P)^(0.2854*(1.0-0.28*q1/100.0))
     ept1=pot1*exp(((3376./tk1)-2.54)*q1/100.0*(1.0+0.81*q1/100.0)) ;
     ept1@description = "0se"
     ept1@units = "K"  
     copy_VarCoords(tc_plane,ept1) ; assign coordinates
;*计算水汽通量散度
  es=(6.112*exp(17.67*(tk-273.15)/(tk-29.65)))
     q=rh_plane*(0.62197*es/(P-es))/100.
   wvflux= q*1000*spd/9.87 A" l- N/ f+ Z% f! m7 Z
     wvflux@description = "moisture flux"

     copy_VarCoords(tc_plane,wvflux) ; assign coordinatesMeteorologi
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   do it = 27,27,1             ; TIME LOOPwww.mnmuc.org9 q4 z4 E7 d4 \7 ?
; do it = 0,ntimes-1,2             ; TIME LOOP
     print("Working on time: " + times(it) )
     res@TimeLabel = times(it)   ;

    do level = 3,3                 ; LOOP OVER LEVELS
       do level = 0,nlevels-1                 ; LOOP OVER LEVELS

      pressure = pressure_levels(level)
       ; Add some level info to the plot
     res@PlotLevelID = pressure + " hPa"
                 
         opts = res                          
      opts@cnLineColor = "Red"
         opts@ContourParameters = (/ 2.0 /)
        opts@cnInfoLabelOrthogonalPosF = 0.07  ; offset second label information
      opts@gsnContourLineThicknessesScale = 2.0
         contour_theatx = wrf_contour(a,wks,ept1(it,level,:,undefined,opts)
   contour_tc = wrf_contour(a,wks,tc_plane(it,level,:,undefined,opts)
         delete(opts)
       ; Plotting options for wvflux         
         opts = res                          
      opts@cnFillOn = True  
        opts@pmLabelBarOrthogonalPosF = -0.1
         contour_wvflux= wrf_contour(a,wks,wvflux(it,level,:,:),opts)

       ; Plotting options for Wind Speed                ( r3 r% c% `* g8 p
         opts = res                          $ c7 r' t; }: v$ A: F
         opts@cnLineColor = "MediumSeaGreen"
         opts@ContourParameters = (/ 10. /)
         opts@cnInfoLabelOrthogonalPosF = 0.07  ;
         opts@gsnContourLineThicknessesScale = 3.0
         contour_spd = wrf_contour(a,wks,spd(it,level,:,:),opts)
         delete(opts)

         opts = res         
   opts@FieldTitle = "Wind"   ; overwrite Field TitleMeteorological Numerical Model Union of China (MNMUC) 中国气象数值模式联盟4 y3 J3 Y. Z% n5 i' Y9 w: N
         opts@NumVectors = 47       ; wind barb density
     vector = wrf_vector(a,wks,u_plane(it,level,:,:),v_plane(it,level,:,:),opts)
       delete(opts)

         opts_z = res                          
         opts_z@cnLineColor = "Blue"
     opts_z@gsnContourLineThicknessesScale = 3.0

         opts_z@ContourParameters = (/ 20.0 /)
         contour_height = wrf_contour(a,wks,z_plane(it,level,:,:),opts_z)
         plot = wrf_map_overlays(a,wks,(/contour_wvflux,contour_theatx,contour_height , \
                                     vector/),pltres,mpres)
   
         delete(opts_z)


end do        ; END OF TIME LOOP

freekiller

实在不好意思，我是外行人！但保证log(x)中的x大于0就要以273为界限啊，为什么要以263为界限呢？

你是不是外行不是你的错，你说我就是你的不对了

     a1 = where(tk .gt. 263.0, 0.622*6.11*exp(17.26*(tk-273.16)/(tk-35.86)), \
                     0.622*6.11*exp(21.87*(tk-273.16)/(tk-7.66)))   ;

这句话的意思就是如果tk>263K, 计算水面饱和水汽，tk<263 K时，计算冰面饱和水汽压。


至于为什么是263 K （-10℃）,这个根据个人的研究对象不同而不同，有的人取-20℃， 如cloudsat数据反演，有的取37℃，有的取-40℃

freekiller

e1需要取对数，log(x)中的x必须大于0，防止x小于0导致计算出问题
翻翻大气物理，你就明白了

juliekxkl

freekiller 发表于 2013-1-3 20:04

                               
登录/注册后可看大图

e1需要取对数，log(x)中的x必须大于0，防止x小于0导致计算出问题
翻翻大气物理，你就明白了

实在不好意思，我是外行人！但保证log(x)中的x大于0就要以273为界限啊，为什么要以263为界限呢？grads里面的'define es=6.112*exp(17.67*(TMPprs-273.15)/(TMPprs-29.65))'
   'define qs=0.62197*es/(prs-0.378*es)'
    'define q=qs*RHprs/100'
   'define e=prs*q/(0.62197+q)+(1e-10)'这几行就相当于
    a1 = where(tk .gt. 263.0, 0.622*6.11*exp(17.26*(tk-273.16)/(tk-35.86)), \
                      0.622*6.11*exp(21.87*(tk-273.16)/(tk-7.66)))   ;
      b1= where(tk .gt. 263.0, P-0.278*exp(17.26*(tk-273.16)/(tk-35.86)),\
                  P-0.278*exp(21.87*(tk-273.16)/(tk-7.66)))   ;
     qs1=a1/b1
q1=qs1*rh_plane
      e1=P*q1/100./(0.62197+q1/100.0)
当tk>263时选where的第二项，但tk-35.86和grads里面的tk-29.65又不一致了，是不是不对啊？要是tk<263的时候选择的where中的第三项，这个公式在哪里可以找到？还请高人指示啊

juliekxkl

freekiller 发表于 2013-1-4 13:54

                               
登录/注册后可看大图

实在不好意思，我是外行人！但保证log(x)中的x大于0就要以273为界限啊，为什么要以263为界限呢？

你是不 ...

太谢谢您的指点了！可是“你说我就是你的不对了”？我没有说你什么啊，感谢你还来不及呢

我是吉星

学西学习，谢谢

菜鸟的梦想

{:5_235:}{:5_235:}{:5_235:}{:5_235:}{:5_235:}{:5_235:}{:5_235:}{:5_235:}

weather3699

{:5_213:}{:5_213:}{:5_213:}

weather3699

我试一试呢

virl

学习，谢谢