请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 5298|回复: 1

[其他] 【求助】以前可以运行出结果的湿位涡,现在全是缺测值

[复制链接]

新浪微博达人勋

发表于 2021-11-29 17:01:34 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

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

搜狗截图20211129170212.jpg
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2022-5-8 00:19:12 来自手机 | 显示全部楼层
请问解决了吗?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表