求助湿Q矢量散度有关问题

绿棕

这是平台里面找的计算湿Q矢量散度的程序，但是对照公式感觉好像少点东西，但又不知如何修改，请问有没有老师可以赐教

'reinit'
'sdfopen D:\202509-10\hgt.2025.nc'
'sdfopen D:\202509-10\uwnd.2025.nc'
'sdfopen D:\202509-10\vwnd.2025.nc'

'set lev 1000 500'
'set lon 90 150'
'set lat 10 60'
'set t 270 295'

'define hgtc=hgt.1'
'define u=uwnd.2'
'define v=vwnd.3'

'define u1=ave(u,t=274,t=291,1)'
'define v1=ave(v,t=274,t=291,1)'

*************计算湿Q矢量散度
'define ug= -9.8*cdiff(hgtc,y)/(111177*cdiff(lat,y))'
'define vg= 9.8*cdiff(hgtc,x)/(111177*cdiff(lon,x)*cos(lat*3.1415/180))'
'define dugdx=cdiff(ug,x)/(111177*cdiff(lon,x)*cos(lat*3.1415/180))'
'define dvgdx=cdiff(vg,x)/(111177*cdiff(lon,x)*cos(lat*3.1415/180))'
'define dtdx=cdiff(hgtc,x)/(111177*cdiff(lon,x)*cos(lat*3.1415/180))'
'define dugdy= cdiff(ug,y)/(111177*cdiff(lat,y))'
'define dvgdy= cdiff(vg,y)/(111177*cdiff(lat,y))'
'define dtdy= cdiff(hgtc,y)/(111177*cdiff(lat,y))'
'define qx=(-287.05/lev)*(dugdx*dtdx+dvgdx*dtdy)*100000000000000'
'define qy=(-287.05/lev)*(dugdy*dtdx+dvgdy*dtdy)*100000000000000'
'define div=hdivg(qx,qy)*100000'

'set mpdset cnworld'
'set map 1 1 1'

'set vpage 0 10 0 8'
'set parea 1 9.8 1 8'
'set grads off'
'set grid off'
'set xlopts 1 4.5 0.19'
'set ylopts 1 4.5 0.19'

'set lon 105 130'
'set lat 25 45'

'set gxout shaded'
'D:\gradscase\grey2.gs'
'set cmax 0'
'd ave(div,t=274,t=291,1)'
'cbarn 1 0'

'set gxout barb'
'set ccolor 1'
'd 2.5*u1;2.5*v1'

'printim D:\202509-10\Q700.png white'
'gxprint D:\202509-10\Q700.pdf x1000 y800 white'
'c'