【求助】用fnl资料求整层水汽通量，结果总是出问题，不知道是哪里出错了

瑾珂

; 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)
      do j = 6,10;201701:8-12;201801:6-10;201812:112-116;
        file_list = addfile(FILE(j), "r")
        time=str_get_cols(FILE(j),40,50)
        p = file_list->lv_ISBL0
        lev=file_list->lv_ISBL0
        num_p=num(p)
        lon = file_list->lon_0
        lat = file_list->lat_0
        tmp = file_list->TMP_P0_L100_GLL0(:,:,:)
        u   = file_list->UGRD_P0_L100_GLL0(:,:,:)
        v   = file_list->VGRD_P0_L100_GLL0(:,:,:)
        rh=file_list->RH_P0_L100_GLL0(:,:,:)
        ps = file_list->PRES_P0_L1_GLL0(:,:)
        rh@_FillValue = -999.0
        es = new((/num_p,181,360/),float,rh@_FillValue)
        es=(6.112*exp((17.67*(tmp-273.16))/(tmp-35.86))) ;saturation vapor pressure
        qs = new((/num_p,181,360/),float,rh@_FillValue)
        do i =0,30
        qs(i,:,:)=(0.62197*es(i,:,:)/(p(i)-0.378*es(i,:,:)))   ;saturation specific humidity
        end do
        q=qs*rh/100             ;specific humidity
        qu=1.0/9.8*u*q
        qv=1.0/9.8*v*q
        copy_VarCoords(tmp,qu)
        copy_VarCoords(tmp,qv)
        qu!0         = "lev"
        qu!1         = "lat"
        qu!2         = "lon"
        qu&lon = lon
        qu&lat = lat
        qu&lev = lev
        lon@units = "degrees_east"
        lat@units = "degrees_north"
        lev@units = "pa"

        qv!0         = "lev"
        qv!1         = "lat"
        qv!2         = "lon"
        qv&lon = lon
        qv&lat = lat
        qv&lev = lev
        lon@units = "degrees_east"
        lat@units = "degrees_north"
        lev@units = "pa"
        linlog=1
        psfc   = 1013.
        pbot=1000
        ptop=100      
        g=9.8
        vint_qu=vibeta(p(::-1),qu(lat|:,lon|:,lev|::-1),linlog,ps,pbot,ptop)/g*100  
        copy_VarMeta(qu(0,:,:),vint_qu)
        vint_qv=vibeta(p(::-1),qv(lat|:,lon|:,lev|::-1),linlog,ps,pbot,ptop)/g*100  
        copy_VarMeta(qv(0,:,:),vint_qv)
        mm=sqrt(vint_qu^2+vint_qv^2)*10^(3)
        copy_VarMeta(vint_qu,mm)
        mmmax=max(mm)
        mmmin=min(mm)
        print(mmmax)
        print(mmmin)


      wks =gsn_open_wks("png","L:/DJBY/ca_saf/201801/waterflu/vapour_flux_"+time)
      gsn_define_colormap(wks,"BlAqGrYeOrReVi200")   
      vcres = True
      vcres@gsnAddCyclic  = False
      vcres@gsnDraw      = False                        ; don't draw yet
      vcres@gsnFrame     = False                        ; don't advance frame yet
      vcres@mpDataSetName         = "Earth..3"   ; This new database contains  
      vcres@mpDataBaseVersion     = "MediumRes"  ; Medium resolution database
      vcres@mpOutlineOn           = True         ; Turn on map outlines
      vcres@mpProjection          = "CylindricalEquidistant"
      vcres@mpGeophysicalLineThicknessF = 0.7      ; double the thickness of geophysical boundaries
      vcres@mpNationalLineThicknessF    = 1.0      ; double the thickness of national boundaries
      vcres@pmTickMarkDisplayMode = "Always"
      vcres@mpMinLatF         = 20                        
      vcres@mpMaxLatF         = 30                        
      vcres@mpMinLonF         = 96                        
      vcres@mpMaxLonF         = 108                                                      
      vcres@lbLabelBarOn = True
      vcres@vcRefMagnitudeF         = 10.0              ; make vectors larger
      vcres@vcRefLengthF            = 0.0018            ; ref vec length
      vcres@vcGlyphStyle            = "CurlyVector"    ; turn on curly vectors
      vcres@vcMinDistanceF          = 0.017            ; thin out vectors
      vcres@vcRefAnnoOn = False
      vcres@vcMonoFillArrowFillColor  = True
      vcres@vcLineArrowThicknessF  = 2.0
      vcres@vcRefAnnoOn  =False
      vcres@vcVectorDrawOrder= "PostDraw"
      vcres@tiMainString          ="vaper_flux : kg/(m*s)"
      vcres@gsnRightString ="Rammasun(1409)"
      plot=gsn_csm_vector_map(wks,vint_qu,vint_qv,vcres)  ; create plot








      res                 = True                    ; plot mods desired
      res@gsnDraw         = False                   ; don't draw yet
      res@gsnFrame        = False                   ; don't advance frame yet
      res@cnFillOn        = True                    ; turn on color
      res@gsnSpreadColors = True                    ; use full colormap
      res@gsnSpreadColorStart = 0
      res@gsnSpreadColorEnd   = 199
      res@cnLinesOn       = False                   ; turn off contour lines
      res@cnLineLabelsOn  = False                   ; tuen off line labels
      res@cnInfoLabelOn =False
      res@gsnLeftString=""
      res@gsnRightString=""
      res@cnLevelSelectionMode     = "ManualLevels"   ; manual contour levels
      res@cnMinLevelValF           = 0          ; minimum level
      res@cnMaxLevelValF           = 50         ; maximum level
      res@lbLabelBarOn = True
      res@lbLabelStride       = 2         ; plot every other colar bar label
      res@lbOrientation        = "vertical"         ; vertical label bars
      res@lbLeftMarginF =-0.5
      ;m=smth9(m,0.5,0.25,False)
      plot1 = gsn_csm_contour(wks,mm,res)
      overlay(plot,plot1)
      frame(wks)  


    end do  
end

   

                               
登录/注册后可看大图