荷兰农气模型WOFOST&ORYZA2000逐日气象数据（CABOWE）文件生成代码

topmad

     好久没有发原创的东西了，根本原因是水平有限，百无一能，看到论坛的大牛们都在晒自己的大作，亚历山大啊，只有埋头苦干，奋起直追了，闲话少叙，回归正题。
     今天发的是一个瓦赫宁根派模型（主要是WOFOST 和 ORYZA2000）的天气文件生成代码，由于模型是事先做好的，说句实话，这些模型的用户体验确实不好，有些东西总感觉不那么舒服顺畅，比如里面所用的逐日气象资料，相对来说，格式的要求十分严格，稍微错一点，最后来一个 simulation filed，那个抓狂的心啊...
     先给大家show一下模型自带的天气文件的格式，看着很简单，如果看了user‘s manual  ，就知道这里面的纠结之处了

1. *---------------------------------------------------------------------------*
   
2. * Station name: Swifterbant, Netherlands
   
3. * Year: 1977
   
4. * Author: Peter Uithol -99.000: NIL VALUE
   
5. * Source: Natuur- en Weerkunde (J.G. Krom)
   
6. * Comments: Missing values are exchanged with daily weather
   
7. * data from synoptic station de Bilt, Netherlands
   
8. * Measured vapour pressure at day 130 was 1.19 kPa
   
9. * this was above saturated vapour pressure 1.16 kPa
   
10. * so the measured vapour pressure was substituted
    
11. * with the saturated vapour pressure.
    
12. *
    
13. * Longitude: 5 38 E, latitude: 52 34 N, altitude: -4 m.
    
14. *
    
15. * Column Daily value
    
16. * 1 station number
    
17. * 2 year
    
18. * 3 day
    
19. * 4 irradiation (kJ m-2 d-1)
    
20. * 5 minimum temperature (degrees Celsius)
    
21. * 6 maximum temperature (degrees Celsius)
    
22. * 7 early morning vapour pressure (kPa)
    
23. * 8 mean wind speed (height: 2 m) (m s-1)
    
24. * 9 precipitation (mm d-1)
    
25. *
    
26. ** WCCDESCRIPTION=Netherlands, Swifterbant
    
27. ** WCCFORMAT=2
    
28. ** WCCYEARNR=1977
    
29. *---------------------------------------------------------------------------*
    
30. 5.63 52.57 -4. -0.18 -0.55
    
31. 2 1977 1 1410. 2.0 9.0 0.810 6.1 2.5
    
32. 2 1977 2 2330. -0.1 5.0 0.650 2.6 0.0
    
33. 2 1977 3 810. 0.3 4.7 0.720 4.6 2.7
    
34. 2 1977 4 1820. -1.0 4.4 0.680 1.7 0.9
    
35. 
    
36. .....

复制代码

    上面的就是模型自带的CABOWE文件，如果将该模式本地化，还要根据当地的气象资料，替换这里面的内容，包括温湿度，风速，水汽压，常规的做法是在excel里面按照格式要求将*-------------*下面的内容做好，然后粘贴到文本中，如果像我做的内容，研究n个省里面nn个站的话，而且还是十余年的内容（一个站一年就要生成一个文件啊），那这样的方法显然就不行了，只有死磕的写个程序了，最初以为很简单，按照格式写一下就行了，直到动了手才知道，里面很复杂，平年闰年，缺测值，地理信息，台站信息，为方便日后理解，我还在台站信息里加了中文台站名，后来又证明模式不认识中文，又加了拼音，还有一些非观测内容比如辐射量，水汽压等等，都要在中间计算，慢慢改了好久，才把最终的文件都做了出来，跟大家一起分享，一个人的水平毕竟有限，希望大家多多之处里面的错误和不足（没有删除测试以及说明行，希望便于大家理解，多多批评）。

     首先，进行数据说明，格式如下，一个省份一个文件


最左边是站号，具体每一列代表什么意思，见下图（示例，着重看看要素代码）


由于还需要在最终的作物文件中写入台站的地理信息，因此又把每个省的台站信息进行了统计，格式如图


每个省的台站数量不同，平闰年的天数不同，为了高效，编程过程中使用了动态数组

贴代码了

1. program JR
   
2. implicit none
   
3. integer i,j,line,m,l,line2,k,sl,sp !line为数据文件的行数 line2 台站信息文件的描述,sl为站名的字符串长度，sp为台站名的长度
   
4. character hehe !求取行数所用
   
5. character*5 sitenumber
   
6. character*2 temp !后缀名转换的中间文件
   
7. character*8 filename !引入汉语站名，欲作为文件名，由于wofost不认可，只将其作为文件内部的注释用
   
8. character*8,allocatable::chinesename(:) !引入了汉语站名 为更方便
   
9. integer,allocatable::provsite(:) !省台站数组
   
10. character*3 lastname !后缀名
    
11. !character*5 sitename !原定以为文件名，由于wofost特殊要求，暂弃用
    
12. character*20 factsite !实际的站名 为汉语拼音，非站号
    
13. real corr !纬度的度数与弧度之间的转换系数
    
14. real lati
    
15. real,allocatable::dat(:,:) !数据文件数组
    
16. real,allocatable::geodata(:,:) ! 台站地理信息 第一列为纬度，第二列精度，第三列海拔高度
    
17. integer,allocatable::site(:) !台站号数组
    
18. character*12,allocatable::actname(:) !台站名数组
    
19. character*8 prov(4)
    
20. REAL LONGTI,heigh !精度 高度
    
21. integer,allocatable::jday(:) !日叙述
    
22. integer*4::result
    
23. integer,external::systemqq
    
24. integer :: time(0:1),timeOff
    
25. real,allocatable::dr(:),detla(:),ws(:),ra(:),rs(:),nn(:),rss(:) !通过日照时数计算辐射的中间变量，ws不是windspeed，需注意
    
26. corr=3.14/180 !公式所用纬度为弧度，corr为换算系数
    
27. CALL SYSTEM_CLOCk(TIME(0))
    
28. data prov/'Hebei','Jiangsu','Henan','Shandong'/
    
29. result=systemqq("mkdir D:\SURFDATA")
    
30. !data provsite/58027,58040,58138,58141,58144,58150,58238,58241,58251,58259,58265,58343,58345,58358/
    
31. do sp=1,4
    
32. result=systemqq("mkdir D:\SURFDATA"//trim(prov(sp))//"")
    
33. line2=0
    
34. open(10,file=''//trim(prov(sp))//'data.txt',status='old')
    
35. do while(10)
    
36. read(10,*,end=200) hehe
    
37. line2=line2+1
    
38. enddo
    
39. 200 continue
    
40. allocate(geodata(line2,3))
    
41. allocate(site(line2))
    
42. allocate(actname(line2))
    
43. allocate(provsite(line2))
    
44. allocate(chinesename(line2))
    
45. rewind(10)
    
46. do i=1,line2
    
47. read(10,*) site(i),actname(i),(geodata(i,j),j=1,3),chinesename(i)
    
48. end do
    
49. close(10)
    
50. do m=1,line2
    
51. provsite(m)=site(m)
    
52. enddo
    
53. ! 测试行
    
54. !do i=1,line2
    
55. ! print*, provsite(i),site(i),geodata(i,1),geodata(i,3),actname(i)
    
56. !enddo
    
57. open (11,file=''//trim(prov(sp))//'.txt',status='old')
    
58. line=0
    
59. do while(11)
    
60. read(11,*,end=100) hehe
    
61. line=line+1
    
62. end do
    
63. 100 continue
    
64. allocate(dat(line-1,11))
    
65. allocate(jday(line-1))
    
66. allocate(dr(line-1))
    
67. allocate(detla(line-1))
    
68. allocate(ws(line-1))
    
69. allocate(ra(line-1))
    
70. allocate(rs(line-1))
    
71. allocate(nn(line-1))
    
72. allocate(rss(line-1))
    
73. !allocate(longti(line))
    
74. !allocate(heigh(line))
    
75. !allocate(factsite(line))
    
76. !allocate(lati(line))
    
77. ! print*,'haha' 测试行
    
78. rewind(11)
    
79. read(11,*) hehe
    
80. do i=1,line-1
    
81. read(11,*) (dat(i,j),j=1,11)
    
82. end do
    
83. close(11)
    
84. 
    
85. !@@@@@@@@@@@@@@@@@@@@@日序计算@@@@@@@@@@@@@@@@@@
    
86. do i=1,line-1
    
87. jday(i)=int(275*(dat(i,3)/9)-30+dat(i,4))-2
    
88. if((mod(int(dat(i,2)),4)==0.and.mod(int(dat(i,2)),100)/=0.or.mod(int(dat(i,2)),400)==0))then
    
89. if(dat(i,3)>2) then
    
90. jday(i)=jday(i)+1
    
91. else
    
92. jday(i)=jday(i)+2
    
93. endif
    
94. elseif(dat(i,3)<3) then
    
95. jday(i)=jday(i)+2
    
96. end if
    
97. enddo
    
98. !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
    
99. 
    
100. do i=1,line-1
     
101. !@@@@@@@@@@@@@不同台站不同纬度，纬度筛选@@@@@@@
     
102. do m=1,line2
     
103. if (int(dat(i,1))==site(m)) then
     
104. lati=geodata(m,1)*corr
     
105. else
     
106. endif
     
107. enddo
     
108. !日照时数数据检测，若缺测，用线性内插法订正
     
109. if(dat(i,11)>=32700) then
     
110. dat(i,11)=(dat(i-1,11)+dat(i+1,11))/2
     
111. endif
     
112. 
     
113. !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
     
114. !@@@@@@@@@@@@@@@计算每日的辐射@@@@@@@@@@@@@@@@
     
115. detla(i)=0.4209*sin(jday(i)*0.0172-1.39)
     
116. dr(i)=1+0.033*cos(0.0172*jday(i))
     
117. ws(i)=acos(-1*tan(lati)*tan(detla(i)))
     
118. ra(i)=37.6*dr(i)*((ws(i))*sin(lati)*sin(detla(i))+cos(lati)*cos(detla(i))*sin(ws(i)))
     
119. nn(i)=24*ws(i)/3.14
     
120. rss(i)=(0.25+(0.5)*dat(i,11)/(10*nn(i)))*ra(i)
     
121. !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
     
122. !print*,rss(i)
     
123. enddo !this do is for do i=1,line
     
124. !stop
     
125. 
     
126. !@@@@@@@@@@@@@@@@@@创建wofost天气文件@@@@@@@@@@@@@@@@@@
     
127. 
     
128. !@@@@@@@@数据订正@@@@@@@
     
129. do i=1,line-1 ! 降水缺测数据筛选，若缺测或痕量，定义为0
     
130. if (dat(i,9)>=32000) then
     
131. dat(i,9)=0
     
132. endif
     
133. enddo
     
134. 
     
135. !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
     
136. do m=1,line2
     
137. do k=1,line2
     
138. if (provsite(m)==site(k)) then !根据台站号进行地理信息的赋值
     
139. lati=geodata(k,1)
     
140. longti=geodata(k,2)
     
141. heigh=geodata(k,3)
     
142. factsite=actname(k)
     
143. filename=chinesename(k)
     
144. write(sitenumber,'(I5)') site(k)
     
145. else
     
146. endif
     
147. enddo
     
148. do l=2000,2011
     
149. if (l-2000<10) then
     
150. write(temp,'(I1)') l-2000
     
151. lastname='00'//trim(temp)//''
     
152. else
     
153. write(temp,'(i2)') l-2000
     
154. lastname='0'//trim(temp)//''
     
155. endif
     
156. open(int(m*l/2),file='D:\surfdata\'//TRIM(PROV(SP))//'\'//trim(factsite)//'1.'//trim(lastname)//'',status='replace')
     
157. sl=len(trim(factsite))
     
158. print*,sl, factsite
     
159. 
     
160. 
     
161. !@@@@@@@@@@@@@@@@@@@ WOFOST头文件创建@@@@@@@@@@@@@@@@@@@@@@
     
162. 
     
163. 
     
164. write(int(m*l/2),"(A63)") '*-------------------------------------------------------------*'
     
165. write(int(m*l/2),"(A15,A8,2X,A5,a17)") '* Station name:',trim(filename),SITENUMBER,','//trim(prov(sp))//'-China'
     
166. WRITE(int(m*l/2),"(A7,I5)") '* Year:',l
     
167. write(int(m*l/2),"(A16,A47)") '* Author: TOPMAD','-99.000: NIL VALUE'
     
168. WRITE(int(m*l/2),"(a15)") '* Source: NUIST'
     
169. WRITE(int(m*l/2),"(A63)") '* Comments: missing data is cauculated with interpoltion method'
     
170. WRITE(int(m*l/2),"(a1)") '*'
     
171. write(int(m*l/2),"(A54)") '* 缺测的日照时数通过前后两日的日照时数之和的平均值订正'
     
172. write(int(m*l/2),"(A31)") '* 痕量以及缺测降水数据均订正为0'
     
173. WRITE(int(m*l/2),"(a1)") '*'
     
174. WRITE(int(m*l/2),"(a)") '* visit bbs.06climate.com get more information aboust us'
     
175. WRITE(int(m*l/2),"(A12,I3,1X,I2,A14,I2,1X,I2,1X,A11,F7.1,A5)") '* Longitude:',INT(longti),INT((LONGTI-INT(LONGTI))*60),'E, latitude: ',INT(lati),INT((LATI-INT(LATI))*60),'N, altitude:',heigh,'m.'
     
176. WRITE(int(m*l/2),"(a1)") '*'
     
177. write(int(m*l/2),"(A21)") '* Column Daily value'
     
178. write(int(m*l/2),"(A24)") '* 1 station number'
     
179. write(int(m*l/2),"(A14)") '* 2 year'
     
180. write(int(m*l/2),"(A13)") '* 3 day'
     
181. write(int(m*l/2),"(A52)") '* 4 irradiation (kJ m-2 d-1)'
     
182. write(int(m*l/2),"(A57)") '* 5 minimum temperature (degrees Celsius)'
     
183. write(int(m*l/2),"(A57)") '* 6 maximum temperature (degrees Celsius)'
     
184. write(int(m*l/2),"(A45)") '* 7 early morning vapour pressure (kPa)'
     
185. write(int(m*l/2),"(A47)") '* 8 mean wind speed (height: 2 m) (m s-1)'
     
186. write(int(m*l/2),"(A48)") '* 9 precipitation (mm d-1)'
     
187. WRITE(int(m*l/2),"(a1)") '*'
     
188. ! if (sl==4.and.sp==1) then
     
189. write(int(m*l/2),"(A18,A,A12)") '** WCCDESCRIPTION=',trim(factsite),','//trim(prov(sp))//'-China'
     
190. WRITE(int(m*l/2),"(a14)") '** WCCFORMAT=2'
     
191. write(int(m*l/2),"(A13,I4)") '** WCCYEARNR=',l
     
192. write(int(m*l/2),"(A63)") '*-------------------------------------------------------------*'
     
193. write(int(m*l/2),"(f7.2,2x,f7.2,2x,f7.2,2x,f7.2,1x,f5.1)") LONGTI,LATI,HEIGH,-0.25,-0.5
     
194. 
     
195. 
     
196. ! wofoso可以用日照时数或者辐射，如果用日照时数，上面的0.25 0.5就要是正的，如果用辐射，就要是0或负数
     
197. do i=1,line-1
     
198. if((int(dat(i,1))==provsite(m)).and.(int(dat(i,2))==l))then
     
199. write(int(m*l/2),'(A1,2X,I4,1X,I3,1X,F10.2,1X,F6.2,1X,F6.2,1X,F6.2,1X,F5.2,1X,F6.2)') '1',INT(dat(i,2)),jday(i),rss(i)*1000,DAT(I,7)/10,DAT(i,6)/10,dat(i,8)/100,(dat(i,10)/10)*4.87/log(678-5.42),dat(i,9)/10
     
200. endif
     
201. enddo !this enddo is for i=1,line
     
202. close(int(m*l/2))
     
203. end do !for m=1,line
     
204. end do ! for l=2000,2011
     
205. !@@@@@@@@@@@释放@@@@@@@@@@@
     
206. deallocate(dat)
     
207. deallocate(jday)
     
208. deallocate(dr)
     
209. deallocate(detla)
     
210. deallocate(ws)
     
211. deallocate(ra)
     
212. deallocate(rs)
     
213. deallocate(nn)
     
214. deallocate(rss)
     
215. deallocate(provsite)
     
216. deallocate(geodata)
     
217. deallocate(site)
     
218. deallocate(actname)
     
219. deallocate(chinesename)
     
220. enddo !this end do is for do sp=1,4(省份循环)
     
221. !@@@@@@@@@用时统计@@@@@@@@@@@
     
222. CALL SYSTEM_CLOCk(TIME(1))
     
223. TIMEOFF=TIME(1)-TIME(0)
     
224. write(*,'(A15,1X,I3,A1)') 'TIME CONSUME：',timeoff,'S'
     
225. end
     

复制代码

这就是完成4个省的cabowe文件的代码，仔细看其实不难，就是相对繁琐一些，如果需要，将里面的代码稍加改动，便可变为己用

希望大家多对立面的任何错误和不足进行批评指正，老五在此先有礼了

zsmnantong

也太长了吧

topmad

zsmnantong 发表于 2011-12-4 10:16

                               
登录/注册后可看大图

也太长了吧

写的水平不行 你看光定义变量跟说明占了多少行

棉花糖

topmad 发表于 2011-12-4 10:17

                               
登录/注册后可看大图

写的水平不行 你看光定义变量跟说明占了多少行

厉害多多学习

mofangbao

老五俨然已经成熟不少了啊

qxtlyf

zhoudog 发表于 2011-12-4 10:13

                               
登录/注册后可看大图

大侠 厉害~~~ 千万别谦虚~~ 我们这些菜鸟还要仰仗您的大作呢

太牛了，学学

topmad

mofangbao 发表于 2011-12-4 11:33

                               
登录/注册后可看大图

老五俨然已经成熟不少了啊

清风大侠指点的到位啊..

时光飞逝

谢谢分享！！！！！！！！！！！！！！！！！1

冬日情愫+

楼主太厉害了！

婧女其姝rat

好东西！mark一下
话说这个WOFOST好用不？靠谱吗