利用GrADs中的aave和mag函数求200公里范围的垂直风切变的问题

pengyulong

1. *本程序用于提取1000hpa-300hpa的垂直分切边
   
2. 'reinit'
   
3. 'open /mnt/usb1/Typoon_Project/fnl00_14.ctl'
   
4. 'set mproj off'
   
5. latstr=readfile('latstr.txt')
   
6. lonstr=readfile('lonstr.txt')
   
7. tstr=readfile('dateNO.txt')
   
8. m=1
   
9. while(m<9075);*num of latstr
   
10. nowlat=subwrd(latstr,m)
    
11. nowlon=subwrd(lonstr,m)
    
12. nowt=subwrd(tstr,m)
    
13. 'set lat ' nowlat
    
14. 'set lon ' nowlon
    
15. 'set t ' nowt
    
16. *-----------------1000hpa-300hpa 200km内的垂直风切变---------------------
    
17. 'set z 1'
    
18. 'display UGRDprs'
    
19.   line=sublin(result,1)
    
20.   u1000=subwrd(line,4)
    
21. 'display VGRDprs'
    
22.   line=sublin(result,1)
    
23.   v1000=subwrd(line,4)
    
24. 'set z 16'
    
25. 'display UGRDprs'
    
26.   line=sublin(result,1)
    
27.   u300=subwrd(line,4)
    
28. 'display VGRDprs'
    
29.   line=sublin(result,1)
    
30.   v300=subwrd(line,4)
    
31. 'define du='u1000-u300''
    
32. 'define dv='v1000-v300''
    
33. 'display aave(mag(du,dv),lon='nowlon-2',lon='nowlon+2',lat='nowlat-2',lat='nowlat+2')'
    
34.   line=sublin(result,1)
    
35.   shr1=subwrd(line,4)
    
36.   write("/mnt/usb1/Typoon_Project/OriginData/shr1_200.txt",line,append)
    
37.   if (rc=0); break; endif
    
38. m=m+1
    
39. endwhile
    
40. function readfile(FileName)
    
41. Lines=''
    
42. while(1);
    
43.     r=read(FileName)
    
44.     temp=subwrd(r,1);*第一行仅有一个单词,为状态码,读取无错误应状态码为0
    
45.     if (temp != 0);break;endif;*如果读取有误或读到EOF就退出
    
46.     line=sublin(r,2);*提取该行字符串
    
47.     Lines=Lines' ' line ' '
    
48. endwhile
    
49. *返回文本字符串
    
50. return(Lines)
    
51. ;
    
52. 'quit'

*本程序用于提取1000hpa-300hpa的垂直分切边<br /> 'reinit'<br /> 'open /mnt/usb1/Typoon_Project/fnl00_14.ctl'<br /> 'set mproj off'<br /> latstr=readfile('latstr.txt')<br /> lonstr=readfile('lonstr.txt')<br /> tstr=readfile('dateNO.txt')<br /> m=1<br /> while(m<9075);*num of latstr<br /> nowlat=subwrd(latstr,m)<br /> nowlon=subwrd(lonstr,m)<br /> nowt=subwrd(tstr,m)<br /> 'set lat ' nowlat<br /> 'set lon ' nowlon<br /> 'set t ' nowt<br /> *-----------------1000hpa-300hpa 200km内的垂直风切变---------------------<br /> 'set z 1'<br /> 'display UGRDprs'<br />   line=sublin(result,1)<br />   u1000=subwrd(line,4)<br /> 'display VGRDprs'<br />   line=sublin(result,1)<br />   v1000=subwrd(line,4)<br /> 'set z 16'<br /> 'display UGRDprs'<br />   line=sublin(result,1)<br />   u300=subwrd(line,4)<br /> 'display VGRDprs'<br />   line=sublin(result,1)<br />   v300=subwrd(line,4)<br /> 'define du='u1000-u300''<br /> 'define dv='v1000-v300''<br /> 'display aave(mag(du,dv),lon='nowlon-2',lon='nowlon+2',lat='nowlat-2',lat='nowlat+2')'<br />   line=sublin(result,1)<br />   shr1=subwrd(line,4)<br />   write("/mnt/usb1/Typoon_Project/OriginData/shr1_200.txt",line,append)<br />   if (rc=0); break; endif<br /> m=m+1<br /> endwhile<br /> function readfile(FileName)<br /> Lines=''<br /> while(1);<br />     r=read(FileName)<br />     temp=subwrd(r,1);*第一行仅有一个单词,为状态码,读取无错误应状态码为0<br />     if (temp != 0);break;endif;*如果读取有误或读到EOF就退出<br />     line=sublin(r,2);*提取该行字符串<br />     Lines=Lines' ' line ' '<br /> endwhile<br /> *返回文本字符串<br /> return(Lines)<br /> ;<br /> 'quit'

  代码中三个文本文件"dateNO.txt"、"latstr.txt","lonstr.txt"都是9074行的数值文件。分别代表日期，纬度，经度。代码都能正常运行，但问题是输出的文件中的数值均是9.8*10^8.这表示根本没有按照display aave(mag(du,dv),lon='nowlon-2',lon='nowlon+2',lat='nowlat-2',lat='nowlat+2')'去求，后面调试程序发现的错误出现在mag(du,dv)上面，但是我却不知道如何修改？望各位大侠帮忙！！

pengyulong

问题已经解决，应该先设置垂直风切变的公式，再设置经纬度和时间。失误啊。

sunflower

楼主给力，弱弱地问一下，请问您用的垂直风切变公式是什么样的呢

pengyulong

垂直风切变，我用这个公式。
'define du = UGRDprs(lev = 850)-UGRDprs(lev = 500)'
'define dv = VGRDprs(lev = 850)-VGRDprs(lev = 500)'
'define shr = mag(du,dv)'*850hpa-500hpa之间的垂直风切变。

pengyulong

sunflower 发表于 2016-1-20 10:27
楼主给力，弱弱地问一下，请问您用的垂直风切变公式是什么样的呢

垂直风切变，我用这个公式。
'define du = UGRDprs(lev = 850)-UGRDprs(lev = 500)'
'define dv = VGRDprs(lev = 850)-VGRDprs(lev = 500)'
'define shr = mag(du,dv)'*850hpa-500hpa之间的垂直风切变。

sunflower

pengyulong 发表于 2016-1-20 15:27
垂直风切变，我用这个公式。
'define du = UGRDprs(lev = 850)-UGRDprs(lev = 500)'
'define dv = VGRD ...

非常感谢楼楼