关于求某段时间变量之和的问题

黄小仙儿

循环部分我用红色标记了，错误提示维度不匹配，我想应该是rh是四维，而rhs是一维，想请问怎么把rsh定义成相对应的四维数组
我试了试直接
do ...
rsh=rh
rsh=rsh+rh
end do
但是转念一想，第一次循环岂不是加了两次rh,实在是很费解
脚本如下：
想求每三个时间序列的变量之和的垂直剖面图，请问大家有没有好的实现方法
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin
;
; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.
  a = addfile("/home/Huanglei/data/d03"+".nc","r")


; We generate plots, but what kind do we prefer?
  type = "pdf"
; type = "pdf"
; type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"plt_CrossSectionsum")
  gsn_define_colormap(wks,"WhBlGrYeRe")

; Set some basic resources
  res = True
  res@MainTitle = "REAL-TIME WRF"

  pltres = True


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  FirstTime = True
  times  = wrf_user_getvar(a,"times",-1) ; get times in the file
  ntimes = dimsizes(times)         ; number of times in the file

  mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file
  nd = dimsizes(mdims)

;---------------------------------------------------------------

  do it = 0,ntimes-1,3             ; TIME LOOP

    print("Working on time: " + times(it) )
    res@TimeLabel = times(it)   ; Set Valid time to use on plots

    tc  = wrf_user_getvar(a,"tc",it)     ; T in C
    rh = wrf_user_getvar(a,"QICE",it)      ; relative humidity
    z   = wrf_user_getvar(a, "z",it)     ; grid point height

    rh = rh*1000.
      rh@units = "g/kg"  

    if ( FirstTime ) then                ; get height info for labels
      zmin = 0.
      zmax = max(z)/1000.
      nz   = floattoint(zmax/2 + 1)
      FirstTime = False
    end if

;---------------------------------------------------------------

    do ip = 1, 3              ; we are doing 3 plots, specifying different start and end points

        opts = True            ; setting start and end times
        plane = new(4,float)

        if(ip .eq. 1) then
          plane = (/  1,98,  199,98  /) ; start x;y & end x;y point
        end if
        if(ip .eq. 2) then
          plane = (/  2,71, 199,71  /) ; start x;y & end x;y point
        end if
        if(ip .eq. 3) then
          plane = (/   1,1, 199,199  /) ; start x;y & end x;y point
        end if
                rsh=0.
                do itt=it,it+2
                rsh=rsh+rh
                end do
        rh_plane = wrf_user_intrp3d(rhs,z,"v",plane,0.,opts)
        tc_plane = wrf_user_intrp3d(tc,z,"v",plane,0.,opts)

        dim = dimsizes(rh_plane)                      ; Find the data span - for use in labels
        zspan = dim(0)


      ; Options for XY Plots
        opts_xy                         = res
        opts_xy@tiYAxisString           = "Height (km)"
        opts_xy@AspectRatio             = 0.75
        opts_xy@cnMissingValPerimOn     = True
        opts_xy@cnMissingValFillColor   = 0
        opts_xy@cnMissingValFillPattern = 11
        opts_xy@tmYLMode                = "Explicit"
        opts_xy@tmYLValues              = fspan(0,zspan,nz)                    ; Create tick marks
        opts_xy@tmYLLabels              = sprintf("%.1f",fspan(zmin,zmax,nz))  ; Create labels
        opts_xy@tiXAxisFontHeightF      = 0.020
        opts_xy@tiYAxisFontHeightF      = 0.020
        opts_xy@tmXBMajorLengthF        = 0.02
        opts_xy@tmYLMajorLengthF        = 0.02
        opts_xy@tmYLLabelFontHeightF    = 0.015
        opts_xy@PlotOrientation         = tc_plane@Orientation


      ; Plotting options for RH
        opts_rh = opts_xy
        opts_rh@pmLabelBarOrthogonalPosF = -0.07
        opts_rh@ContourParameters       = (/ 0.02, 0.24, 0.02 /)
        opts_rh@cnFillOn                = True
      ;  opts_rh@cnFillColors            = (/"White","White","White", \
       ;                                     "White","Chartreuse","Green", \
      ;                                      "Green3","Green4", \
       ;                                     "ForestGreen","PaleGreen4"/)

      ; Plotting options for Temperature
        opts_tc = opts_xy
        opts_tc@cnInfoLabelOrthogonalPosF = 0.00
        opts_tc@ContourParameters  = (/ 5. /)


      ; Get the contour info for the rh and temp
        contour_tc = wrf_contour(a,wks,tc_plane,opts_tc)
        contour_rh = wrf_contour(a,wks,rh_plane,opts_rh)


      ; MAKE PLOTS         
        plot = wrf_overlays(a,wks,(/contour_rh,contour_tc/),pltres)

      ; Delete options and fields, so we don't have carry over
        delete(opts_tc)
        delete(opts_rh)
        delete(tc_plane)
        delete(rh_plane)

    end do  ; make next cross section

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  end do        ; END OF TIME LOOP

end

longlivehj

我稍微改了一下，你看看吧！
我这边wrfout里面没有QICE变量，换用rh试了一下，可以出图！

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin
;
; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.
  a = addfile("/home/Huanglei/data/d03"+".nc","r")


; We generate plots, but what kind do we prefer?
  type = "pdf"
; type = "pdf"
; type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"plt_CrossSectionsum")
  gsn_define_colormap(wks,"WhBlGrYeRe")

; Set some basic resources
  res = True
  res@MainTitle = "REAL-TIME WRF"

  pltres = True


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  FirstTime = True
  times  = wrf_user_getvar(a,"times",-1) ; get times in the file
  ntimes = dimsizes(times)         ; number of times in the file

  mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file
  nd = dimsizes(mdims)

;---------------------------------------------------------------

  z = wrf_user_getvar(a, "z", 0)     ; grid point height

  do it = 0,ntimes-1,3             ; TIME LOOP

    print("Working on time: " + times(it) )
    res@TimeLabel = times(it)   ; Set Valid time to use on plots

    print(it + (/0, 1, 2/))

    tc  = wrf_user_getvar(a,"tc", it + (/0, 1, 2/))     ; T in C
    rh = wrf_user_getvar(a,"QICE", it + (/0, 1, 2/))      ; relative humidity

    rh = rh*1000.
    rh@units = "g/kg"

    tcs = dim_sum_n_Wrap(tc, 0)
    rhs = dim_sum_n_Wrap(rh, 0)

    if ( FirstTime ) then                ; get height info for labels
      zmin = 0.
      zmax = max(z)/1000.
      nz   = floattoint(zmax/2 + 1)
      FirstTime = False
    end if

;---------------------------------------------------------------

    do ip = 1, 3              ; we are doing 3 plots, specifying different start and end points

        opts = True            ; setting start and end times
        plane = new(4,float)

        if(ip .eq. 1) then
          plane = (/  1,98,  199,98  /) ; start x;y & end x;y point
        end if
        if(ip .eq. 2) then
          plane = (/  2,71, 199,71  /) ; start x;y & end x;y point
        end if
        if(ip .eq. 3) then
          plane = (/   1,1, 199,199  /) ; start x;y & end x;y point
        end if

        rh_plane = wrf_user_intrp3d(rhs,z,"v",plane,0.,opts)
        tc_plane = wrf_user_intrp3d(tcs,z,"v",plane,0.,opts)

        dim = dimsizes(rh_plane)                      ; Find the data span - for use in labels
        zspan = dim(0)


      ; Options for XY Plots
        opts_xy                         = res
        opts_xy@tiYAxisString           = "Height (km)"
        opts_xy@AspectRatio             = 0.75
        opts_xy@cnMissingValPerimOn     = True
        opts_xy@cnMissingValFillColor   = 0
        opts_xy@cnMissingValFillPattern = 11
        opts_xy@tmYLMode                = "Explicit"
        opts_xy@tmYLValues              = fspan(0,zspan,nz)                    ; Create tick marks
        opts_xy@tmYLLabels              = sprintf("%.1f",fspan(zmin,zmax,nz))  ; Create labels
        opts_xy@tiXAxisFontHeightF      = 0.020
        opts_xy@tiYAxisFontHeightF      = 0.020
        opts_xy@tmXBMajorLengthF        = 0.02
        opts_xy@tmYLMajorLengthF        = 0.02
        opts_xy@tmYLLabelFontHeightF    = 0.015
        opts_xy@PlotOrientation         = tc_plane@Orientation


      ; Plotting options for RH
        opts_rh = opts_xy
        opts_rh@pmLabelBarOrthogonalPosF = -0.07
        opts_rh@ContourParameters       = (/ 0.02, 0.24, 0.02 /)
        opts_rh@cnFillOn                = True
      ;  opts_rh@cnFillColors            = (/"White","White","White", \
       ;                                     "White","Chartreuse","Green", \
      ;                                      "Green3","Green4", \
       ;                                     "ForestGreen","PaleGreen4"/)

      ; Plotting options for Temperature
        opts_tc = opts_xy
        opts_tc@cnInfoLabelOrthogonalPosF = 0.00
        opts_tc@ContourParameters  = (/ 5. /)


      ; Get the contour info for the rh and temp
        contour_tc = wrf_contour(a,wks,tc_plane,opts_tc)
        contour_rh = wrf_contour(a,wks,rh_plane,opts_rh)


      ; MAKE PLOTS         
        plot = wrf_overlays(a,wks,(/contour_rh,contour_tc/),pltres)

      ; Delete options and fields, so we don't have carry over
        delete(opts_tc)
        delete(opts_rh)
        delete(tc_plane)
        delete(rh_plane)

    end do  ; make next cross section

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  end do        ; END OF TIME LOOP

end

黄小仙儿

哇，原来可以在数据读入上下功夫，明白了。但是我运行到一段时间后出现了这样的错误：

(2)     53
(0)     Working on time: 2010-07-25_06:00:00
(0)     54
(1)     55
(2)     56
fatal:Subscript out of range, error in subscript #0
fatal:Execute: Error occurred at or near line 1028 in file $NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl

fatal:Execute: Error occurred at or near line 53 in file /home/Huanglei/wrf_CrossSectionsum.ncl

longlivehj

黄小仙儿 发表于 2014-4-10 12:29
哇，原来可以在数据读入上下功夫，明白了。但是我运行到一段时间后出现了这样的错误：

(2)     53

因为你的ntimes不能被3整除。
改成do it = 0,ntimes-3,3