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

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 199|回复: 2

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

[复制链接]

新浪微博达人勋

发表于 2024-4-17 15:31:59 | 显示全部楼层 |阅读模式
20金钱
本帖最后由 HANSEN 于 2024-4-17 15:44 编辑

各位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


图1 暴雨WRF模拟

图1 暴雨WRF模拟

图2 诊断的各项垂直分布

图2 诊断的各项垂直分布

图3 涡度方程

图3 涡度方程
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2024-4-17 15:34:03 | 显示全部楼层
说明一下,图2中900hPa以下是地形引起的缺失值,恳请路过的大佬帮帮我{:5_275:}
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2024-4-17 15:41:22 | 显示全部楼层
我反复检查了很多次脚本,没有发现有什么问题,眼睛已经要瞎了
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

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

本版积分规则

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

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

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