WRF进行涡度方程诊断的残差项过大问题

HANSEN

各位UU，最近在研究用WRF输出的高分辨数据进行涡度方程的诊断，一般来说因为摩擦或者次网格误差因素，方程的左右两边通常不等，会有残差项，但是查询一些文献，都指出残差项通常是小项，可以忽略不计。但是我在诊断的过程中发现，残差项的量级几乎与其他几项量级相当了，这让我百思不得其解。
我这里是随便拿了一个暴雨个例做测试，模拟2022年7月15日陕西一次暴雨事件，模式用FNL驱动，只设置一层网格，分辨率12km，垂直方向设置了53层，如图1给出了模拟的暴雨分布，WRF的模拟是成功的。模式时间步长为60s，同时我也输出了最高分辨率60s的wrfout数据，因为考虑到 偏R/偏t（R为相对涡度）希望t尽量趋近于0，我把wrfout的数据首先插值到P坐标系，然后插值成经纬度网格（0.1°），选取了暴雨中心区域计算涡度各项，如图2得到了各项的垂直剖面图，绿色虚线为残差项，紫色局地变化项，红色平流项，橙色散度项，另外还有一些小项不做说明。
问题来了：为什么这个残差项会这么大呢？
正常来说，我分辨率设置的已经非常高了，涡度方程的右边各项应该尽可能等于左边叭，我不能接受残差项的跟其他各项几乎相当！！
是不是我的脚本写错了，跪求各位大佬帮我看看！或者有没有做过这种分析的大佬传授一下经验

脚本如下：

begin

    path = "~/WRF4.2-6/wrf/run/wrfout_d01_2022-07-14_18:00:00"
    a = addfile(path,"r")
    times  = wrf_user_getvar(a,"times",-1) ; get times in the file
    lat2d = a->XLAT(0,:,:)
    lon2d = a->XLONG(0,:,:)
    minlat = min(lat2d)
    maxlat = max(lat2d)
    minlon = min(lon2d)
    maxlon = max(lon2d)

    lon           = fspan(minlon, maxlon, toint((maxlon-minlon)/0.1)+1)
    lon!0         = "lon"
    lon&lon       = lon
    lon@units     = "degrees_east"
    lon@long_name = "Longitude"

    lat           = fspan(minlat, maxlat, toint((maxlat-minlat)/0.1)+1)
    lat!0         = "lat"
    lat&lat       = lat
    lat@units     = "degrees_north"
    lat@long_name = "Latitude"   

    pp = new(13, "float")
    pp = ispan(800, 500, 25)
    pp@units = "hPa"

    duration = 3
    k= 0

    do hh = 599,601

        ua  = wrf_user_getvar(a,"ua",hh)
        va  = wrf_user_getvar(a,"va",hh)
        wa  = wrf_user_getvar(a,"wa",hh)
        prs  = wrf_user_getvar(a,"pressure",hh)

        u1 = wrf_user_interp_level(ua,prs,pp,False)
        u2 = rcm2rgrid_Wrap(lat2d, lon2d, u1, lat, lon, 1)

        delete(prs)
        pp1 = pp * 100.
        prs = conform_dims(dimsizes(u2), pp1, 0)
        printVarSummary(prs)
        print(prs(5,:,:))

        v1 = wrf_user_interp_level(va,prs,pp,False)
        v2 = rcm2rgrid_Wrap(lat2d, lon2d, v1, lat, lon, 1)
        w1 = wrf_user_interp_level(wa,prs,pp,False)
        w2 = rcm2rgrid_Wrap(lat2d, lon2d, w1, lat, lon, 1)

        if (hh.eq.599) then

            dims = dimsizes(u2)
            VOR = new((/duration,dims(0),dims(1),dims(2)/),float)
            UU = new((/duration,dims(0),dims(1),dims(2)/),float)
            VV = new((/duration,dims(0),dims(1),dims(2)/),float)
            WW = new((/duration,dims(0),dims(1),dims(2)/),float)

        end if

        UU(k,:,:,:) = u2
        VV(k,:,:,:) = v2   
        WW(k,:,:,:) = w2

        k = k + 1

    end do
    delete(prs)
    DIV = uv2dv_cfd(UU, VV, lat, lon, 2)
    VOR = uv2vr_cfd(UU, VV, lat, lon, 2)

    dt = 60
    dRdt = center_finite_diff_n(VOR, dt, False, 0, 0)  * 1000000000.

    pp1 = pp * 100.
    prs = conform_dims(dimsizes(u2), pp1, 0)

    nlat = dimsizes(lat)
    nlon = dimsizes(lon)
    nlevel = dimsizes(pp1)
    dx = new((/nlevel,nlat,nlon/), float)
    dy = dx
    f  = dx

    dlon = (lon(1)-lon(0)) * 0.0174533
    do nx = 0,nlat-1  

        dy(:,nx,:) = 6378388.*lat(nx)*0.0174533
        f(:,nx,:) = coriolis_param(lat(nx))

        do ny = 0,nlon-1
            dx(:,nx,ny) = 6378388.*cos(0.0174533*lat(nx))*dlon * ny
        end do

    end do

    delete([/ua,va,wa,u1,u2,v1,v2,w1,w2/])

    do j = 1,1

        div = DIV(j,:,:,:)
        vor = VOR(j,:,:,:)
        u = UU(j,:,:,:)
        v = VV(j,:,:,:)
        w = WW(j,:,:,:)

        vor!0 = "level"
        vor!1 = "lat"
        vor!2 = "lon"
        vor&level = pp
        vor&lat = lat  
        vor&lon = lon


        ;--------- Term 1: -V·▽R ---------
        dRdx = center_finite_diff_n(vor, dx, False, 0, 2)
        dRdy = center_finite_diff_n(vor, dy, False, 0, 1)
        term1 = -(u*dRdx + v*dRdy) * 1000000000.


        ;--------- Term 2: -w*dRdp ---------
        term2 = - w * center_finite_diff_n(vor, prs, False, 0, 0) * 1000000000.


        ;--------- Term 3: -(R+f)*D ---------
        term3 = - (vor+f) * div * 1000000000.


        ;--------- Term 4: -v*dfdy ---------
        dfdy = center_finite_diff_n(f, dy, False, 0, 1)
        term4 = -v*dfdy * 1000000000.

        ;--------- Term 5: dwdy*dudp-dwdx*dvdp ---------
        dwdx = center_finite_diff_n(w, dx, False, 0, 2)
        dwdy = center_finite_diff_n(w, dy, False, 0, 1)
        dudp = center_finite_diff_n(u, prs, False, 0, 0)
        dvdp = center_finite_diff_n(v, prs, False, 0, 0)
        term5 = (dwdy*dudp-dwdx*dvdp) * 1000000000.

        copy_VarCoords(vor, term1)
        copy_VarCoords(vor, term2)
        copy_VarCoords(vor, term3)
        copy_VarCoords(vor, term4)
        copy_VarCoords(vor, term5)
        copy_VarCoords(vor, dRdt(0,:,:,:))

        maxlat1 = 36.
        minlat1 = 35.
        maxlon1 = 110.
        minlon1 = 109.

        dRdt_avg = dim_avg_n_Wrap(dRdt(j,:,{minlat1:maxlat1},{minlon1:maxlon1}), (/1,2/))
        term1_avg = dim_avg_n_Wrap(term1(:,{minlat1:maxlat1},{minlon1:maxlon1}), (/1,2/))
        term2_avg = dim_avg_n_Wrap(term2(:,{minlat1:maxlat1},{minlon1:maxlon1}), (/1,2/))
        term3_avg = dim_avg_n_Wrap(term3(:,{minlat1:maxlat1},{minlon1:maxlon1}), (/1,2/))
        term4_avg = dim_avg_n_Wrap(term4(:,{minlat1:maxlat1},{minlon1:maxlon1}), (/1,2/))
        term5_avg = dim_avg_n_Wrap(term5(:,{minlat1:maxlat1},{minlon1:maxlon1}), (/1,2/))

        resi = dRdt_avg - term1_avg - term2_avg - term3_avg - term4_avg - term5_avg
        total = term1_avg + term2_avg + term3_avg + term4_avg + term5_avg

        copy_VarCoords(vor, resi)

        wks = gsn_open_wks("png", "vertical_vor_eqt_test")
        gsn_define_colormap(wks, "Cat12")

        lev = toint(pp1 / 100.)
        lev1 =lev(::-1)

        res                        = True
        res@gsnFrame               = False
        res@gsnDraw                = False
        res@gsnMaximize            = True

        res@tiMainString           = ""
        res@tiMainFontHeightF      = 0.03
        res@tiXAxisString          = "10~S~-9~N~ "
        res@tiYAxisString          = "Pressure Level (hPa)"
        res@tiXAxisFontHeightF     = 0.015
        res@tiYAxisFontHeightF     = 0.015

        res@trXMinF           = -15
        res@trXMaxF           = 15
        res@trYMaxF           = 800            
        res@trYMinF           = 500         

        res@xyLineThicknesses = 5
        res@xyDashPatterns    = 0

        res@tmYLMode          = "Explicit"
        res@tmYLValues        = lev(0::2)
        res@tmYLLabels        = tostring(lev1(0::2))
        res@tmYLLabelFontHeightF   = 0.015
        res@tmYLMajorLengthF  = 0.015
        res@tmYLMajorOutwardLengthF = 0.015
        res@tmXBLabelFontHeightF   = 0.015


        res@gsnXRefLine = 0
        res@gsnXRefLineDashPattern = 1

        res@tmXTOn = False
        res@tmYROn = False

        res@xyLineColors      = 11
        plot1 = gsn_csm_xy(wks,dRdt_avg,lev(::-1),res)

        res@xyLineColors      = 13
        plot2 = gsn_csm_xy(wks,term1_avg,lev(::-1),res)
        overlay(plot1, plot2)

        res@xyLineColors      = 12
        plot3 = gsn_csm_xy(wks,term2_avg,lev(::-1),res)
        overlay(plot1, plot3)

        res@xyLineColors      = 3
        plot4 = gsn_csm_xy(wks,term3_avg,lev(::-1),res)
        overlay(plot1, plot4)

        res@xyLineColors      = 5
        plot6 = gsn_csm_xy(wks,term4_avg,lev(::-1),res)
        overlay(plot1, plot6)

        res@xyLineColors      = 9
        plot7 = gsn_csm_xy(wks,term5_avg,lev(::-1),res)
        overlay(plot1, plot7)

        res@xyDashPatterns    = 1
        res@xyLineColors      = 7
        plot5 = gsn_csm_xy(wks,resi,lev(::-1),res)
        overlay(plot1, plot5)

        draw(plot1)
        frame(wks)

    end do

end

HANSEN

说明一下，图2中900hPa以下是地形引起的缺失值，恳请路过的大佬帮帮我{:5_275:}

HANSEN

我反复检查了很多次脚本，没有发现有什么问题，眼睛已经要瞎了