正压原始方程组（浅水波方程）fortran源代码

ningking1986

1. !i is the direction of y, j is the direction of x.
   
2. !the firt mesh point is origin of coordinate
   
3. program barotropic
   
4. implicit none
   
5. real(8),parameter :: g=9.8,d=3.e5,dt=1.0,s=0.5
   
6. integer,parameter :: mi=16,mj=20,mk=48
   
7. real(8) z0(mi,mj),z(mi,mj,0:mk),u(mi,mj,0:mk),v(mi,mj,0:mk),rm(mi,mj),f(mi,mj)
   
8. !----------------------------------------------- the value of z,u,v,ku,kv,kz at the time dt/2 is z1,u1,v1,ku1,kv1,kz1
   
9. real(8) z1(mi,mj),u1(mi,mj),v1(mi,mj),kz(mi,mj,0:mk),ku(mi,mj,0:mk),kv(mi,mj,0:mk),kz1(mi,mj),ku1(mi,mj),kv1(mi,mj)
   
10. integer i,j,k
    
11. !----------------------------------------------read the geopotential height field
    
12. open(1,file='data.txt')
    
13. read(1,*)((z0(i,j),j=1,mj),i=1,mi)
    
14. close(1)
    
15. 
    
16. call mf(d,rm,f,mi,mj) !calculate m,f
    
17. forall(i=1:mi,j=1:mj) z0(i,j)=z0(i,j)*10
    
18. !-------------------------------------------------smooth the geopotential filed
    
19. do i=2,mi-1
    
20.   do j=2,mj-1
    
21.     z0(i,j)=z0(i,j)+s/4*(z0(i+1,j)+z0(i-1,j)+z0(i,j+1)+z0(i,j-1)-4*z0(i,j))
    
22.   enddo
    
23. enddo
    
24. !-------------------------------------------------calculate the initial value of u,v,z
    
25. do i=2,mi-1
    
26.   do j=2,mj-1
    
27.     z(i,j,0)=z0(i,j)
    
28.         u(i,j,0)=-rm(i,j)*g/f(i,j)*(z0(i-1,j)-z0(i+1,j))/2/d
    
29.     v(i,j,0)=rm(i,j)*g/f(i,j)*(z0(i,j+1)-z0(i,j-1))/2/d
    
30.   enddo
    
31. enddo
    
32. !-------------------------------------------------calculate the boundary value of u,v,z at the points on the direction of y
    
33. 
    
34. do k=0,mk
    
35.   do i=2,mi-1
    
36.      z(i,1,k)=z0(i,1)
    
37.      z(i,mj,k)=z0(i,mj)
    
38.      u(i,1,k)=-rm(i,1)*g/f(i,1)*(z0(i-1,1)-z0(i+1,1))/2/d      !central difference
    
39.      u(i,mj,k)=-rm(i,mj)*g/f(i,mj)*(z0(i-1,mj)-z0(i+1,mj))/2/d !central difference
    
40.      v(i,1,k)=rm(i,1)*g/f(i,1)*(z0(i,2)-z0(i,1))/d             !forward difference
    
41.      v(i,mj,k)=rm(i,mj)*g/f(i,mj)*(z0(i,mj)-z0(i,mj-1))/d      !backward difference
    
42.   enddo
    
43. !-------------------------------------------------calculate the boundary value of u,v,z at the points on the direction of x
    
44.   do j=2,mj-1
    
45.       z(1,j,k)=z0(1,j)
    
46.       z(mi,j,k)=z0(mi,j)
    
47.       u(1,j,k)=-rm(1,j)*g/f(1,j)*(z0(1,j)-z0(2,j))/d            !backward difference
    
48.       u(mi,j,k)=-rm(mi,j)*g/f(mi,j)*(z0(mi-1,j)-z0(mi,j))/d     !forward difference
    
49.       v(1,j,k)=rm(1,j)*g/f(1,j)*(z0(1,j+1)-z0(1,j-1))/2/d       !central difference
    
50.       v(mi,j,k)=rm(mi,j)*g/f(mi,j)*(z0(mi,j+1)-z0(mi,j-1))/2/d  !central difference
    
51.   enddo
    
52. !-------------------------------------------------calculate the boundary value of u,v,z at the points on the four corners
    
53.   z(1,1,k)=z0(1,1)
    
54.   z(1,mj,k)=z0(1,mj)
    
55.   z(mi,1,k)=z0(mi,1)
    
56.   z(mi,mj,k)=z0(mi,mj)
    
57. !   u(1,1,k)=0.5*(u(2,1,k)+u(1,2,k))
    
58. !   u(1,mj,k)=0.5*(u(2,mj,k)+u(1,mj-1,k))
    
59. !   u(mi,1,k)=0.5*(u(mi-1,1,k)+u(mi,2,k))
    
60. !   u(mi,mj,k)=0.5*(u(mi-1,mj,k)+u(mi,mj-1,k))
    
61. !   v(1,1,k)=0.5*(v(2,1,k)+v(1,2,k))
    
62. !   v(1,mj,k)=0.5*(v(2,mj,k)+v(1,mj-1,k))
    
63. !   v(mi,1,k)=0.5*(v(mi-1,1,k)+v(mi,2,k))
    
64. !   v(mi,mj,k)=0.5*(v(mi-1,mj,k)+v(mi,mj-1,k))
    
65.    u(1,1,k)=-rm(1,1)*g/f(1,1)*(z0(1,1)-z0(2,1))/d               !backward difference
    
66.    u(1,mj,k)=-rm(1,mj)*g/f(1,mj)*(z0(1,mj)-z0(2,mj))/d          !backward difference
    
67.    u(mi,1,k)=-rm(mi,1)*g/f(mi,1)*(z0(mi-1,1)-z0(mi,1))/d          !forward difference
    
68.    u(mi,mj,k)=-rm(mi,mj)*g/f(mi,mj)*(z0(mi-1,mj)-z0(mi,mj))/d     !forward difference
    
69. 
    
70.    v(1,1,k)=rm(1,1)*g/f(1,1)*(z0(1,2)-z0(1,1))/d                !forward difference
    
71.    v(1,mj,k)=rm(1,mj)*g/f(1,mj)*(z0(1,mj)-z0(1,mj-1))/d           !backward difference
    
72.    v(mi,1,k)=rm(mi,1)*g/f(mi,1)*(z0(mi,2)-z0(mi,1))/d           !forward difference
    
73.    v(mi,mj,k)=rm(mi,mj)*g/f(mi,mj)*(z0(mi,mj)-z0(mi,mj-1))/d      !backward difference
    
74. enddo
    
75. !-------------------------------------------------------forward difference and central difference ,the step is half-step
    
76. k=0
    
77. call time_qiancha(u,v,z,u1,v1,z1,rm,f,k,mi,mj,mk)
    
78. call time_zhongyangcha_one(u,v,z,u1,v1,z1,rm,f,k,mi,mj,mk)
    
79. !------------------------------------------------------central difference, the step is one-step
    
80. do k=1,47
    
81.   call time_zhongyangcha(u,v,z,rm,f,k,mi,mj,mk)
    
82.     if(mod(k+1,12)==0)then !temporal smoothing once every 12 hours
    
83.       do i=2,mi-1 !temporal smoothing
    
84.         do j=2,mj-1
    
85.          u(i,j,k)=u(i,j,k)+s/2*(u(i,j,k+1)+u(i,j,k-1)-2*u(i,j,k))
    
86.          v(i,j,k)=v(i,j,k)+s/2*(v(i,j,k+1)+v(i,j,k-1)-2*v(i,j,k))
    
87.          z(i,j,k)=z(i,j,k)+s/2*(z(i,j,k+1)+z(i,j,k-1)-2*z(i,j,k))
    
88.         enddo
    
89.       enddo
    
90.       do i=2,mi-1 !spatial smoothing
    
91.         do j=2,mj-1
    
92.          u(i,j,k)=u(i,j,k)+s/4*(u(i+1,j,k)+u(i-1,j,k)+u(i,j+1,k)+u(i,j-1,k)-4*u(i,j,k))
    
93.          v(i,j,k)=v(i,j,k)+s/4*(v(i+1,j,k)+v(i-1,j,k)+v(i,j+1,k)+v(i,j-1,k)-4*v(i,j,k))
    
94.          z(i,j,k)=z(i,j,k)+s/4*(z(i+1,j,k)+z(i-1,j,k)+z(i,j+1,k)+z(i,j-1,k)-4*z(i,j,k))
    
95.         enddo
    
96.       enddo
    
97.     endif
    
98. enddo
    
99. 
    
100. open(2,file='result.txt')
     
101. open(3,file='10.txt')
     
102. open(4,file='mf.txt')
     
103. do k=0,mk
     
104.   write(2,100)k,'the value of z'
     
105.     do i=1,mi
     
106.       write(2,200)(z(i,j,k),j=1,mj)
     
107.     enddo
     
108. 
     
109.     write(2,100)k,'the value of u'
     
110. 
     
111.     do i=1,mi
     
112.       write(2,200)(u(i,j,k),j=1,mj)
     
113.     enddo
     
114. 
     
115.     write(2,100)k,'the value of v'
     
116. 
     
117.     do i=1,mi
     
118.       write(2,200)(v(i,j,k),j=1,mj)
     
119.     enddo
     
120.        
     
121.         write(3,200)(u(7,j,k),j=1,mj)       
     
122. enddo
     
123.    write(4,*)'the value of m'
     
124.    write(4,200)((rm(i,j),j=1,mj),i=1,mi)
     
125.    write(4,*)'the value of f'
     
126.    write(4,200)((f(i,j),j=1,mj),i=1,mi)
     
127. close(2)
     
128. close(3)
     
129. close(4)
     
130. 100 format(i4,1x,a)
     
131. 200 format(<mj>f12.6)
     
132. end
     
133. ! ----------------------------------------------------------------define the subroutines
     
134. subroutine time_qiancha(u,v,z,u1,v1,z1,rm,f,k,mi,mj,mk)
     
135. real(8),parameter :: g=9.8,d=3.e5,dt=1.0,s=0.5
     
136. real(8) z(mi,mj,0:mk),u(mi,mj,0:mk),v(mi,mj,0:mk),rm(mi,mj),f(mi,mj),kz(mi,mj,0:mk),ku(mi,mj,0:mk),kv(mi,mj,0:mk),z1(mi,mj),u1(mi,mj),v1(mi,mj)
     
137. integer k
     
138. call space(u,v,z,rm,f,k,ku,kv,kz,mi,mj,mk)
     
139. do i=2,mi-1
     
140.   do j=2,mj-1
     
141.     u1(i,j)=u(i,j,k)+dt/2*ku(i,j,k)
     
142.     v1(i,j)=v(i,j,k)+dt/2*kv(i,j,k)
     
143.     z1(i,j)=z(i,j,k)+dt/2*kz(i,j,k)
     
144.   enddo
     
145. enddo
     
146. end subroutine
     
147. !-----------------------------------------------------
     
148. subroutine time_zhongyangcha_one(u,v,z,u1,v1,z1,rm,f,k,mi,mj,mk)
     
149. real(8),parameter :: g=9.8,d=3.e5,dt=1.0,s=0.5
     
150. real(8) z(mi,mj,0:mk),u(mi,mj,0:mk),v(mi,mj,0:mk),z1(mi,mj),u1(mi,mj),v1(mi,mj),rm(mi,mj),f(mi,mj),kz1(mi,mj),ku1(mi,mj),kv1(mi,mj)
     
151. integer k
     
152. call space_one(u1,v1,z1,ku1,kv1,kz1,rm,f,k,mi,mj)
     
153. do i=2,mi-1
     
154.   do j=2,mj-1
     
155.     u(i,j,k+1)=u(i,j,k)+dt*ku1(i,j)
     
156.     v(i,j,k+1)=v(i,j,k)+dt*kv1(i,j)
     
157.     z(i,j,k+1)=z(i,j,k)+dt*kz1(i,j)
     
158.   enddo
     
159. enddo
     
160. end subroutine
     
161. !-------------------------------------------------------
     
162. subroutine time_zhongyangcha(u,v,z,rm,f,k,mi,mj,mk)
     
163. real(8),parameter :: g=9.8,d=3.e5,dt=1.0,s=0.5
     
164. real(8) z(mi,mj,0:mk),u(mi,mj,0:mk),v(mi,mj,0:mk),rm(mi,mj),f(mi,mj),kz(mi,mj,0:mk),ku(mi,mj,0:mk),kv(mi,mj,0:mk)
     
165. integer k
     
166. call space(u,v,z,rm,f,k,ku,kv,kz,mi,mj,mk)
     
167. do i=2,mi-1
     
168.   do j=2,mj-1
     
169.     u(i,j,k+1)=u(i,j,k-1)+dt*ku(i,j,k)
     
170.     v(i,j,k+1)=v(i,j,k-1)+dt*kv(i,j,k)
     
171.     z(i,j,k+1)=z(i,j,k-1)+dt*kz(i,j,k)
     
172.   enddo
     
173. enddo
     
174. end subroutine
     
175. !----------------------------------------------------- the rate of change to u,v and z respectively is ku,kv and kz
     
176. subroutine space(u,v,z,rm,f,k,ku,kv,kz,mi,mj,mk)
     
177. real(8),parameter :: g=9.8,d=3.e5,dt=1.0,s=0.5
     
178. real(8) z(mi,mj,0:mk),u(mi,mj,0:mk),v(mi,mj,0:mk),rm(mi,mj),f(mi,mj),kz(mi,mj,0:mk),ku(mi,mj,0:mk),kv(mi,mj,0:mk),f0
     
179. integer k
     
180. do i=2,mi-1
     
181.   do j=2,mj-1
     
182.     f0=f(i,j)+1.0/2/d*u(i,j,k)*( rm(i-1,j)-rm(i+1,j) )-1.0/2/d*v(i,j,k)*( rm(i,j+1)-rm(i,j-1) )
     
183.     ku(i,j,k)=-rm(i,j)*( (1.0/4/d)*( ( u(i,j+1,k)+u(i,j,k) )*( u(i,j+1,k)-u(i,j,k) )+( u(i,j,k)+u(i,j-1,k) )*( u(i,j,k)-u(i,j-1,k) ) )+(1.0/4/d)*( ( v(i,j,k)+v(i+1,j,k) )*( u(i,j,k)-u(i+1,j,k) )+( v(i-1,j,k)+v(i,j,k) )*( u(i-1,j,k)-u(i,j,k) ) )+g/2/d*( z(i,j+1,k)-z(i,j-1,k) ) )+f0*v(i,j,k)
     
184.     kv(i,j,k)=-rm(i,j)*( (1.0/4/d)*( ( u(i,j+1,k)+u(i,j,k) )*( v(i,j+1,k)-v(i,j,k) )+( u(i,j,k)+u(i,j-1,k) )*( v(i,j,k)-v(i,j-1,k) ) )+(1.0/4/d)*( ( v(i,j,k)+v(i+1,j,k) )*( v(i,j,k)-v(i+1,j,k) )+( v(i-1,j,k)+v(i,j,k) )*( v(i-1,j,k)-v(i,j,k) ) )+g/2/d*( z(i-1,j,k)-z(i+1,j,k) ) )-f0*u(i,j,k)
     
185.     kz(i,j,k)=-rm(i,j)**2*( (1.0/4/d)*( (u(i,j+1,k)+u(i,j,k))*(z(i,j+1,k)/rm(i,j+1)-z(i,j,k)/rm(i,j))+(u(i,j,k)+u(i,j-1,k))*(z(i,j,k)/rm(i,j)-z(i,j-1,k)/rm(i,j-1)) )+(1.0/4/d)*( (v(i-1,j,k)+v(i,j,k))*(z(i-1,j,k)/rm(i-1,j)-z(i,j,k)/rm(i,j))+(v(i,j,k)+v(i+1,j,k))*(z(i,j,k)/rm(i,j)-z(i+1,j,k)/rm(i+1,j)) )+z(i,j,k)/rm(i,j)*( (1.0/2/d)*(u(i,j+1,k)-u(i,j-1,k))+(1.0/2/d)*(v(i-1,j,k)-v(i+1,j,k)) ) )
     
186. !        write(*,*)ku(i,j,k)
     
187.   enddo
     
188. enddo
     
189. end subroutine
     
190. 
     
191. !----------------------------------------------------
     
192. subroutine space_one(u1,v1,z1,ku1,kv1,kz1,rm,f,k,mi,mj)
     
193. real(8),parameter :: g=9.8,d=3.e5,dt=1.0,s=0.5
     
194. real(8) z1(mi,mj),u1(mi,mj),v1(mi,mj),rm(mi,mj),f(mi,mj),kz1(mi,mj),ku1(mi,mj),kv1(mi,mj),f0
     
195. integer k
     
196. do i=2,mi-1
     
197.   do j=2,mj-1
     
198. 
     
199.     ku1(i,j)=-rm(i,j)*( (1.0/4/d)*( ( u1(i,j+1)+u1(i,j) )*( u1(i,j+1)-u1(i,j) )+( u1(i,j)+u1(i,j-1) )*( u1(i,j)-u1(i,j-1) ) )+(1.0/4/d)*( ( v1(i,j)+v1(i+1,j) )*( u1(i-1,j)-u1(i,j) )+( v1(i-1,j)+v1(i,j) )*( u1(i,j)-u1(i+1,j) ) )+g/2/d*( z1(i,j+1)-z1(i,j-1) )+f(i,j)*v1(i,j)+u1(i,j)*(1.0/2/d)*( rm(i-1,j)-rm(i+1,j) )*v1(i,j)-v1(i,j)*(1.0/2/d)*( rm(i,j+1)-rm(i,j-1) )*v1(i,j) )
     
200.         kv1(i,j)=-rm(i,j)*( (1.0/4/d)*( ( u1(i,j+1)+u1(i,j) )*( v1(i,j+1)-v1(i,j) )+( u1(i,j)+u1(i,j-1) )*( v1(i,j)-v1(i,j-1) ) )+(1.0/4/d)*( ( v1(i,j)+v1(i+1,j) )*( v1(i-1,j)-v1(i,j) )+( v1(i-1,j)+v1(i,j) )*( v1(i,j)-v1(i+1,j) ) )+g/2/d*( z1(i-1,j)-z1(i+1,j) )-f(i,j)*u1(i,j)-u1(i,j)*(1.0/2/d)*( rm(i-1,j)-rm(i+1,j) )*u1(i,j)+v1(i,j)*(1.0/2/d)*( rm(i,j+1)-rm(i,j-1) )*u1(i,j) )
     
201.     kz1(i,j)=-rm(i,j)**2*( (1.0/4/d)*( (u1(i,j+1)+u1(i,j))*(z1(i,j+1)/rm(i,j+1)-z1(i,j)/rm(i,j))+(u1(i,j)+u1(i,j-1))*(z1(i,j)/rm(i,j)-z1(i,j-1)/rm(i,j-1)) )+(1.0/4/d)*( (v1(i-1,j)+v1(i,j))*(z1(i-1,j)/rm(i-1,j)-z1(i,j)/rm(i,j))+(v1(i,j)+v1(i+1,j))*(z1(i,j)/rm(i,j)-z1(i+1,j)/rm(i+1,j)) )+z1(i,j)/rm(i,j)*( (1.0/2/d)*(u1(i,j+1)-u1(i,j-1))+(1.0/2/d)*(v1(i-1,j)-v1(i+1,j)) ) )
     
202. 
     
203. !  write(*,*)kz1(i,j)
     
204.   enddo
     
205. enddo
     
206. end subroutine
     
207. !---------------------------------------------------calculate m,f
     
208. subroutine mf(d,rm,f,mi,mj)
     
209. real(8),parameter :: pi=3.1415926535,omiga=7.292115e-5,a=6.371004e6,ak=0.7156,le=11142370
     
210. real(8) ll(mi,mj),rm(mi,mj),l0,fai0,sita(mi,mj),fai(mi,mj),f(mi,mj),d
     
211. !real(8) weidu(mi,mj),theta(mi,mj)
     
212. fai0=60*pi/180
     
213. l0=a*sin(pi/6)/ak
     
214. do i=1,mi
     
215.   do j=1,mj
     
216.     ll(i,j)=(((j-1)*d)**2+(l0+(i-1)*d)**2)**0.5
     
217.     sita(i,j)=2*atan((ll(i,j)/l0)**(1/ak)*tan(0.5*pi/6))
     
218.     rm(i,j)=ak*ll(i,j)/(a*sin(sita(i,j)))
     
219.     fai(i,j)=fai0-(i-1)*d/a
     
220.     f(i,j)=2*omiga*sin(fai(i,j))
     
221. !  ll(i,j)=(-4+j-1)**2+(9+i)**2
     
222. !  ll(i,j)=sqrt(ll(i,j))*d
     
223. !   weidu(i,j)=(le**(2/ak)-ll(i,j)**(2/ak))/(le**(2/ak)+ll(i,j)**(2/ak))
     
224. !  
     
225. ! weidu(i,j)=asind(weidu(i,j))
     
226. ! theta(i,j)=90.0-weidu(i,j)
     
227. ! rm(i,j)=ak*ll(i,j)/(a*sind(theta(i,j)))
     
228. ! f(i,j) =2.0*omiga*cosd(theta(i,j))
     
229. !        write(*,*)rm(i,5)
     
230.   enddo
     
231. enddo
     
232. !  open(10,file='rm.txt')
     
233. ! do i=1,mi
     
234. ! write(10,300)(rm(i,j),j=1,mj)
     
235. ! enddo
     
236. ! 300 format(20f8.4)
     
237. end subroutine
     

复制代码

ningking1986

把相应得测试数据也顺便传上来吧，五个贡献有点儿贵啊。

南国的佳人

赞！谢谢分享

冒失鬼

随便看看～

tsl08

谢谢，感谢分享~

huangbi_h

感谢楼主分享，楼主辛苦了