求K指数和比湿的gs代码

vmboring

求K指数的gs代码：
'define t850=air.1(lev=850)-273.15'
'define t700=air.1(lev=700)-273.15'
'define t500=air.1(lev=500)-273.15'
'define td850=t850-(14.55+0.114*t850)*(1-0.01*rhum.2(lev=850))+pow((2.5+0.007*t850)*(1-0.01*rhum.2(lev=850)),3)+(15.9+0.117*t850)*pow((1-0.01*rhum.2(lev=850)),14))'
'define td700=t700-(14.55+0.114*t700)*(1-0.01*rhum.2(lev=700))+pow((2.5+0.007*t700)*(1-0.01*rhum.2(lev=700)),3)+(15.9+0.117*t700)*pow((1-0.01*rhum.2(lev=700)),14))'
'define kindex=t850-t500+td850-(t700-td700)'
据气温和相对湿度求比湿：
'define es=6.1078*exp(17.2693882*(air.1-273.16)/(air.1-35.86))'
'define esi=6.1078*exp(21.8745584*(air.1-276.16)/(air.1-7.66))'
if(air.1>273.15)
'define er=es'
else
'define er=esi'
endif
'define shum=(622*rhum.2*er/100)/(lev-0.378*er)'