两点间垂直剖面新问题

追风的阳光

选取wrf输出结果（nX*nY*nZ））的子区域做计算，得到一个nX1*nY1*nZ1的三维数组，将此数组赋予wrf坐标。用wrf_user_ll_to_ij得到母网格上两点的X，Y值，再分别减去子网格相对于母网格的网格数x，y（x=nX-nX1,y=nY-ny2）.得到相对于子网格的x，y值。以此做剖面，请问这样可行吗？脚本如下，这样做取得的剖面不是我想要的，是这种情况下没法用wrf_user_intrp3d函数吗，不知道哪里出问题了。谢谢大家

;****************************************************************
  diri="./forc_ai/"
  diri2="./ADAB-out/"
  files = systemfunc("ls "+diri+"*_solution.txt")
  files2= systemfunc("ls "+diri2+"*_rhot_base_avg.txt")
  Time=(/0,12,24,36,48,60/)

;*****************get resource for xlon and xlat***********************************************
  f = addfile("./wrfout_d03_2015-06-11_21:00:00","r")
;**********************************************************************************************
diri3="./WRFOUTPUT/"
  filepath = systemfunc("ls "+diri3+"wrfout_d03*")
  print(filepath)
  f1  = addfiles(filepath,"r")

;**********************************************************************************************
wks=gsn_open_wks("ps","cross-section-4terms")
gsn_define_colormap(wks,"MPL_jet")
contour_term=new(6,graphic)

;**********************************************************************************************

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

  xlat = wrf_user_getvar(f, "XLAT",0)
  xlon = wrf_user_getvar(f, "XLONG",0)
  printVarSummary(xlat)
  u    =wrf_user_getvar(f, "ua", 0)
  z    =wrf_user_getvar(f, "z", 0)
  XLAT=xlat(115:371,62:446)
  XLON=xlon(115:371,62:446)
;*****************levs***********************************************************************
zbot=650
dz=50
nlev=257
ztop=zbot+dz*(nlev-1)
levs=fspan(zbot,ztop,nlev)  
z1 =wrf_user_intrp3d(z,z,"h",levs,0.,False)

Z=z1(:,115:371,62:446)

printVarSummary(Z)



;*********************************************************************************************
Term_avg_plane=new((/6,100,151/),float)
;*********************get cross-section definition*****************************
     lon_list=(/-97.75,-97.55,-97.35,-97.25,-97.05,-96.85,-96.65/)
     lat_list=(/38.5,38.7,38.9,39.1,39.3,39.5,39.7/)
     distance=(/0,29,49,67,85,101,130/)
     distance2=(/"-75","-40","-20","0","20","40","75"/)
;**********************************resource***************************************************************************
  res = True
    res@vpWidthF =0.8
    res@gsnDraw=False
    res@gsnFrame=False
    res@lbLabelBarOn        = False
    res@cnInfoLabelOn       = False
    res@lbTitleString ="Distance(km)"
    res@tiXAxisString ="Distance(km)"
    res@tiYAxisString ="height(m)"
    ; res@tiXAxisPosition ="Right"
    ; res@tiXAxisOffsetYF =0.03
    ; res@tiXAxisOffsetXF =0.012
    res@tmYLMode      ="Explicit"
    res@tmXBMode      ="EXplicit"
    res@tmYLValues    =(/0,3,9,15,21,27,33,39,45,51,57,63,69,75,81,87/)
    res@tmYLLabels    =(/"600","1000","2000","3000","4000","5000","6000","7000","8000","9000","10000","11000","12000","13000","14000","15000"/)
    res@tmXBValues    =distance
    res@tmXBLabels    =distance2

    res@tiXAxisFontHeightF      = 0.020
    res@tiYAxisFontHeightF      = 0.020
    res@tmXBMajorLengthF        = 0.02
    res@tmYLMajorLengthF        = 0.02
    res@tmYLLabelFontHeightF    = 0.015
    ;res@PlotOrientation         = tc_plane@Orientation
;*******************************res_tc********************************
    res_term=res
    res_term@cnFillOn                = True
    res_term@cnLinesOn               =False
    res_term@cnLineDashPattern       = 0
    res_term@cnLevelSelectionMode = "ManualLevels"
    res_term@cnMinLevelValF        = -10     ; Min contour
    res_term@cnMaxLevelValF        = 10.0     ; Max contour
    res_term@cnLevelSpacingF       = 0.2
    res_term@cnInfoLabelSide = "Top"
    res_term@cnInfoLabelOn      =False
    res_term@cnInfoLabelPerimOn = False
    res_term@cnInfoLabelOrthogonalPosF = -0.00005
    res_term@cnFillPalette ="MPL_jet"
;*******************************************************************res_vc0=res_vc
    res_term0=res_term
    res_term0@gsnLeftString   = "2000 UTC"
    res_term1=res_term
    res_term1@gsnLeftString   = "2100 UTC"
    res_term2=res_term
    res_term2@gsnLeftString   = "2200 UTC"
    res_term3=res_term
    res_term3@gsnLeftString   = "2300 UTC"
    res_term4=res_term
    res_term4@gsnLeftString   = "0000 UTC"
    res_term5=res_term
    res_term5@gsnLeftString   = "0100 UTC"


;********************************************************************************************
   do k=0,5,1
    Term_sum_plane=new((/100,151/),float)
    Term_sum_plane=0
   it=Time(k)
  ;it=0
  Term=asciiread(files(it), (/257,257,385/), "float")
  Rhot=asciiread(files2(it),(/257/),"float")
  Rhot2 = conform_dims(dimsizes(Term),Rhot,(/0/))
  printVarSummary(Rhot2)
  Term=Term/Rhot2
  copy_VarMeta(u, Term)
; copy_VarMeta(u, Rhot)
  Term!0="lev"
  Term!1="degrees_north"
  Term!2="degrees_east"

  Term@lev=Z(:,0,0)
  Term@lat2d=XLAT
  Term@lon2d=XLON

  ; Rhot&bottom_top=Z(:,0,0)
  ; Rhot&south_north=XLAT(:,0)
  ; Rhot&west_east=XLON(0,:)

  printVarSummary(Term)
; printVarSummary(Rhot)  
;***********************************************************************************************
  angle1 =atan2(110,85)
        print(angle1)
   do n=0,3,1
   print("processing n="+n)



        if (n.eq.0.or.2.or.3)
         minlat   =lat_list(n)-0.7*cos(angle1);0.7005
         maxlat   =lat_list(n)+0.7*cos(angle1)
         minlon   =lon_list(n)-0.7*sin(angle1)
         maxlon   =lon_list(n)+0.7*sin(angle1)
         print(minlat+" "+maxlat+" "+minlon+" "+maxlon)
        end if

        if (n.eq.1)   
               minlat   =lat_list(n)-0.7*cos(angle1);0.7
               maxlat   =lat_list(n)+0.7*cos(angle1)
               minlon   =lon_list(n)-0.7*sin(angle1)
               maxlon   =lon_list(n)+0.7*sin(angle1)
               print(minlat+" "+maxlat+" "+minlon+" "+maxlon)
        end if


         loc=wrf_user_ll_to_ij(f1, (/minlon,maxlon/),(/minlat,maxlat/), True)
         print(loc(:,:))
         plane= new(4,float)
         loc(0,0)=loc(0,0)-115
         loc(0,1)=loc(0,0)-115
         loc(1,0)=loc(1,0)-62
         loc(1,1)=loc(1,1)-62
         plane= (/loc(0,0),loc(1,1),loc(0,1),loc(1,0)/);(/ mdims(nd-1)/2, mdims(nd-2)/2 /)    ; pivot point is center of domain (x,y)



;**********************************************************************************        
         ; opts = True
         ; angle =atan2(loc(0,0)-loc(0,1),loc(1,1)-loc(1,0))
         ;  X_plane = wrf_user_intrp2d(XLON,plane,angle,opts)
         ;  print(X_plane)
         ; X_desc = "Distance(150km)"
;**********************************************************************************

  opts=True

  z_plane   =wrf_user_intrp3d(z,z,"v",plane,0.0,opts)
  Term_plane   =wrf_user_intrp3d(Term,Z,"v",plane,0.0,opts)
printVarSummary(z_plane)
  printVarSummary(Term_plane)

;******************do the sum****************************************************************************************************
  Term_sum_plane=Term_sum_plane+Term_plane

end do
;******************do the average************************************************************************************************

  Term_avg_plane(k,:,:)=Term_sum_plane/4
end do
;**************************************************************************************************************
; ;**************************set vertical layer************************************************************************
         zmin = 650.
         zmax = 15000;6000.                          ; We are only interested in the first 6km
         nz   = floattoint((zmax-zmin)/100 + 1)        


         b = ind(z_plane(:,0) .gt. zmax )
         c = ind(z_plane(:,0).gt.zmin)
           zmax_pos = b(0) - 1
           zmin_pos = c(0) - 1
           if (abs(z_plane(zmax_pos,0)-zmax).lt.abs(z_plane(zmax_pos+1,0)-zmax)) then
             zspan = b(0) - 1
           else
             zspan = b(0)
           end if

           delete(z_plane)
           delete(b)
           FirstTime = False
;***********************************************************************************************************************  

do k=0,5,1
    if k.eq.0
        res_term=res_term0
    end if
    if k.eq.1
        res_term=res_term1
    end if
    if k.eq.2
        res_term=res_term2
    end if
    if k.eq.3
        res_term=res_term3
    end if
    if k.eq.4
        res_term=res_term4
    end if
    if k.eq.5
        res_term=res_term5
    end if
contour_term(k)=gsn_csm_contour(wks, 1.e3*Term_avg_plane(k,zmin_pos:zmax_pos,:), res_term)

;print(Term_avg_plane(0,10:12,:)+" "+Term_avg_plane(2,10:12,:))
  end do


;***********************************************************************************************

;***********************wks**********************************************************************************   

     pres              =True
     pres@gsnPanelRowSpec = True
     pres@gsnPanelLabelBar    = True                ; add common colorbar
     pres@lbLabelFontHeightF  = 0.007               ; make labels smaller
     pres@lbTitleString=""
     pres@lbTitleFontHeightF =0.01
     pres@lbTitlePosition ="Bottom"
     gsn_panel(wks,contour_term,(/3,3/),pres)
     frame(wks)


    end


;