计算相对湿度的fortran问题

zhksoft

今天有个工作要根据相对湿度计算湿球温度的问题，搜了一下论坛，看到了这个贴子

fortran编程 数组越界问题 编译的是湿球温度迭代算法 谁 帮帮我！！！
http://bbs.06climate.com/forum.p ... 8&fromuid=35781
(出处: 气象家园)


调试运行了一下，结果发现运行结果和实际的相差很大，就根据地面规范附2

附录2　湿度参量的计算公式

1、饱和水汽压    在一定温度下，空气中的水汽与相毗连的水或冰平面处于相变平衡时湿空气中的水汽压。    饱和水汽压采用世界气象组织推荐的戈夫-格雷奇（Goff-Gratch）公式。    ⑴ 纯水平液面饱和水汽压的计算公式    logEW=10.79574(1-T1/T)-5.02800log(T/T1)+1.50475×10-4［1-10-8.2969（T/ T1-1）］+0.42873×10-3(104.76955(1-T1/T)-1)］+0.78614       式中，EW：纯水平液面饱和水汽压（hPa）；T1=273.16K（水的三相点温度）；T=273.15+t℃(绝对温度K)。    ⑵ 纯水平冰面饱和水汽压的计算公式    LogEi=-9.09685[T1/(T-1)]-3.56654Log(T1/T)+0.87682[1-  T1/(T-1)]+0.78614    式中，Ei：水平冰面饱和水汽压（hPa）；T1和T同上。    2、水汽压    ⑴ 用干湿球温度求空气中水汽压的计算公式    e=Etw-APh(t-tw)    中，e：水汽压（hPa）；Etw：湿球温度tw所对应的纯水平液面的饱和水汽压，湿球结冰且湿球温度低于0℃时，为纯水平冰面的饱和水汽压；A：干湿表系数（℃-1），由干湿表类型、通风速度及湿球结冰与否而定，其值见干湿表系数表；Ph：本站气压（hPa）；t：干球温度（℃）；tW：湿球温度（℃）。               附表2.1 干湿表系数表                                                                                   干湿表类型及通风速度                                    Ai×10-3(℃-1)                                    湿球未结冰                                    湿球结冰                                    通风干湿表（通风速度2.5m/s）                                    0.662                                    0.584                                    球状干湿表(通风速度0.4m/s)                                    0.857                                    0.756                                    柱状干湿表(通风速度0.4m/s)                                    0.815                                    0.719                                    现用百叶箱球状干湿表(通风速度0.8m/s)                                    0.7947                                    0.7947                   ⑵ 当使用湿敏电容、毛发表或湿度计等直接测得相对湿度时，由相对湿度求水汽压公式    　　　e=U×EW/100   　　式中，U：相对湿度（%）；e：水汽压（hPa）；EW：干球温度t所对应的纯水平液面饱和水汽压（hPa）。    3、相对湿度    ⑴ 使用干湿球温度表测湿时，空气中相对湿度的计算公式    　　　　　U=（E/EW）×100%    式中，U相对湿度（%）；E水汽压（hPa）；Ew干球温度t所对应的纯水平液面（或冰面）饱和水汽压（hPa）。    ⑵ 使用毛发湿度表（计）测湿时，空气中相对湿度的计算公式    　　　　　Y=bo+b1X+b2X2+b3X3    式中，Y：经毛发湿度表（计）订正后的相对湿度（%）；X：毛发湿度表（计）读数（%）； bo、b1、b2、b3：回归多项式系数，即毛发湿度表（计）的订正系数。

想先用正推的办法算相对湿度，结果算出来 的相对湿度 水汽压 露点等都不对，我的代码如下：

1. program mian
   
2. !声明变量
   
3. real::e !水汽压
   
4. real::t !干球温度
   
5. real::p !本站气压
   
6. real,parameter::T1=273.16
   
7. real::t0
   
8. integer,parameter::N=1500
   
9. real::tw(0:n) !湿球温度估算值序列
   
10. real::etw,etw_g,etw_s !饱和水汽压
    
11. !real,parameter::a=8.15*10**(-4)
    
12. real,parameter::A=7.947*10**(-4)
    
13. real::ej(0:n) !水汽压的计算值
    
14. real::err(0:n) !计算误差
    
15. real::ermin !计算误差的最小值
    
16. 
    
17. real::u,t_g,t_s,td
    
18. 
    
19. !输入变量
    
20. write(*,*)"水汽压为"
    
21. !read(*,*)e
    
22. write(*,*)"干球温度为"
    
23. !read(*,*)t
    
24. write(*,*)"本站气压为"
    
25. !read(*,*)p
    
26. e=14.7
    
27. t_g=25.5
    
28. t_s=17.3
    
29. !t_g=26.3
    
30. !t_s=18.3
    
31. !t_g=31.4
    
32. !t_s=19.1
    
33. P=998.6
    
34. 
    
35. t0=273.15+t_s  !湿球
    
36. etw_s=10**(10.79574*(1-T1/t0)-5.02800*log(t0/T1)+1.50475*(10**(-4))*(1-10**((-8.2969*(t0/T1-1))))+0.42873*(10**(-3))*(10**(4.76955*(1-T1/t0))-1)+0.78614)
    
37. t0=273.15+t_g  !干球
    
38. etw_g=10**(10.79574*(1-T1/t0)-5.02800*log(t0/T1)+1.50475*(10**(-4))*(1-10**((-8.2969*(t0/T1-1))))+0.42873*(10**(-3))*(10**(4.76955*(1-T1/t0))-1)+0.78614)
    
39. e=etw_s-A*P*(t_g-t_s)
    
40. td=(243.92*log(etw_g/6.1078))/(7.69-log(etw_g/6.1078))
    
41. 
    
42. write(*,*)"=======\n湿球ETw="
    
43. write(*,*)etw_s
    
44. write(*,*)"=======\n干球ETw="
    
45. write(*,*)etw_g
    
46. write(*,*)"=======\n水汽压e="
    
47. write(*,*)e
    
48. 
    
49. u=(e/etw_g)*100
    
50. write(*,*)"=======\n相对湿度u="
    
51. write(*,*)u
    
52. write(*,*)"=======\n露点温度="
    
53. write(*,*)td
    
54. 
    
55. stop
    
56. end
    

program mian<br /> !声明变量<br /> real::e !水汽压<br /> real::t !干球温度<br /> real::p !本站气压<br /> real,parameter::T1=273.16<br /> real::t0<br /> integer,parameter::N=1500<br /> real::tw(0:n) !湿球温度估算值序列<br /> real::etw,etw_g,etw_s !饱和水汽压<br /> !real,parameter::a=8.15*10**(-4)<br /> real,parameter::A=7.947*10**(-4)<br /> real::ej(0:n) !水汽压的计算值<br /> real::err(0:n) !计算误差<br /> real::ermin !计算误差的最小值<br /> <br /> real::u,t_g,t_s,td<br /> <br /> !输入变量 <br /> write(*,*)"水汽压为"<br /> !read(*,*)e<br /> write(*,*)"干球温度为"<br /> !read(*,*)t<br /> write(*,*)"本站气压为"<br /> !read(*,*)p<br /> e=14.7<br /> t_g=25.5<br /> t_s=17.3<br /> !t_g=26.3<br /> !t_s=18.3<br /> !t_g=31.4<br /> !t_s=19.1<br /> P=998.6<br /> <br /> t0=273.15+t_s  !湿球<br /> etw_s=10**(10.79574*(1-T1/t0)-5.02800*log(t0/T1)+1.50475*(10**(-4))*(1-10**((-8.2969*(t0/T1-1))))+0.42873*(10**(-3))*(10**(4.76955*(1-T1/t0))-1)+0.78614)<br /> t0=273.15+t_g  !干球<br /> etw_g=10**(10.79574*(1-T1/t0)-5.02800*log(t0/T1)+1.50475*(10**(-4))*(1-10**((-8.2969*(t0/T1-1))))+0.42873*(10**(-3))*(10**(4.76955*(1-T1/t0))-1)+0.78614)<br /> e=etw_s-A*P*(t_g-t_s)<br /> td=(243.92*log(etw_g/6.1078))/(7.69-log(etw_g/6.1078))<br /> <br /> write(*,*)"=======\n湿球ETw="<br /> write(*,*)etw_s<br /> write(*,*)"=======\n干球ETw="<br /> write(*,*)etw_g<br /> write(*,*)"=======\n水汽压e="<br /> write(*,*)e<br /> <br /> u=(e/etw_g)*100<br /> write(*,*)"=======\n相对湿度u="<br /> write(*,*)u<br /> write(*,*)"=======\n露点温度="<br /> write(*,*)td<br /> <br /> stop<br /> end<br />

各位大神看看哪里出了问题。

zhksoft

没人关注？自顶一下

zhksoft

怎么回事呢？大家对此不感兴趣？

言深深

这个东西不是大家不感兴趣，而是太枯燥了。
建议自己再检查一下，逐步print出来，看看到底哪一步出错，可能容易一些。
单单从程序上面看，你能跑出来结果，就说明程序是通的。结果不对，只可能是过程算的不对，比如方程是否写错，参数是否带错，运算过程是否注意到IN规则，等等。把一些过程变量print出来看看，问题应该就能够解决了。

sxnsq8888

现在用的少了，不太记得了

zhksoft

费了N多脑细胞，终于搞定了{:5_275:}

爱在深秋

关注一下吧

DAISYQH

请问楼主相对湿度直接*100就行了吗？不用输出百分号？
如果想把相对湿度按照70%这种形式输出来该怎么做呀？

samvelly

关注一下吧