Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

interpol_wind_nests.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 11 years ago
sources for flexwrf v3.1
File size: 9.2 KB

Line
1	!**********************************************************************
2	! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
3	! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
4	! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
5	! *
6	! This file is part of FLEXPART. *
7	! *
8	! FLEXPART is free software: you can redistribute it and/or modify *
9	! it under the terms of the GNU General Public License as published by*
10	! the Free Software Foundation, either version 3 of the License, or *
11	! (at your option) any later version. *
12	! *
13	! FLEXPART is distributed in the hope that it will be useful, *
14	! but WITHOUT ANY WARRANTY; without even the implied warranty of *
15	! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
16	! GNU General Public License for more details. *
17	! *
18	! You should have received a copy of the GNU General Public License *
19	! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
20	!**********************************************************************
21
22	subroutine interpol_wind_nests(itime,xt,yt,zt, &
23	uprof,vprof,wprof, usigprof,vsigprof,wsigprof, &
24	rhoprof,rhogradprof, tkeprof,pttprof, &
25	u,v,w,usig,vsig,wsig,pvi, &
26	p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2, &
27	ix,jy,ixp,jyp,ngrid,indz,indzp, depoindicator, &
28	indzindicator,mu,mv)
29	! i i i i
30	!*****************************************************************************
31	! *
32	! This subroutine interpolates the wind data to current trajectory position.*
33	! *
34	! Author: A. Stohl *
35	! *
36	! 16 December 1997 *
37	! 16 December 1997 *
38	! *
39	! Revision March 2005 by AST : all output variables in common block cal- *
40	! culation of standard deviation done in this *
41	! routine rather than subroutine call in order *
42	! to save computation time *
43	! *
44	!*****************************************************************************
45	! *
46	! Variables: *
47	! u,v,w wind components *
48	! itime [s] current temporal position *
49	! memtime(3) [s] times of the wind fields in memory *
50	! xt,yt,zt coordinates position for which wind data shall be *
51	! calculated *
52	! *
53	! Constants: *
54	! *
55	!*****************************************************************************
56
57	use par_mod
58	use com_mod
59	! use interpol_mod
60
61	implicit none
62
63	integer :: itime
64	real :: xt,yt,zt
65
66	! Auxiliary variables needed for interpolation
67	real :: dz1,dz2,dz
68	real :: u1(2),v1(2),w1(2),uh(2),vh(2),wh(2)
69	real :: usl,vsl,wsl,usq,vsq,wsq,xaux
70	integer :: i,m,n,indexh,indzh
71	real,parameter :: eps=1.0e-30
72
73	real :: uprof(nzmax),vprof(nzmax),wprof(nzmax)
74	real :: usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax)
75	real :: rhoprof(nzmax),rhogradprof(nzmax)
76	real :: tkeprof(nzmax),pttprof(nzmax)
77	real :: u,v,w,usig,vsig,wsig,pvi,mu,mv
78
79	real :: p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2
80	integer :: ix,jy,ixp,jyp,ngrid,indz,indzp
81	logical :: depoindicator(maxspec)
82	logical :: indzindicator(nzmax)
83
84
85	!********************************************
86	! Multilinear interpolation in time and space
87	!********************************************
88
89	! Determine the lower left corner and its distance to the current position
90	!*************************************************************************
91
92	ddx=xt-real(ix)
93	ddy=yt-real(jy)
94	rddx=1.-ddx
95	rddy=1.-ddy
96	p1=rddx*rddy
97	p2=ddx*rddy
98	p3=rddx*ddy
99	p4=ddx*ddy
100
101	mu =p1*m_xn(ix ,jy ,1,ngrid) &
102	+ p2*m_xn(ixp,jy ,1,ngrid) &
103	+ p3*m_xn(ix ,jyp,1,ngrid) &
104	+ p4*m_xn(ixp,jyp,1,ngrid)
105	mv =p1*m_yn(ix ,jy ,1,ngrid) &
106	+ p2*m_yn(ixp,jy ,1,ngrid) &
107	+ p3*m_yn(ix ,jyp,1,ngrid) &
108	+ p4*m_yn(ixp,jyp,1,ngrid)
109
110
111	! Calculate variables for time interpolation
112	!*******************************************
113
114	dt1=real(itime-memtime(1))
115	dt2=real(memtime(2)-itime)
116	dtt=1./(dt1+dt2)
117
118	! Determine the level below the current position for u,v
119	!*******************************************************
120
121	do i=2,nz
122	if (height(i).gt.zt) then
123	indz=i-1
124	goto 6
125	endif
126	end do
127	6 continue
128
129
130	! Vertical distance to the level below and above current position
131	!****************************************************************
132
133	dz=1./(height(indz+1)-height(indz))
134	dz1=(zt-height(indz))*dz
135	dz2=(height(indz+1)-zt)*dz
136
137
138	!**********************************************************************
139	! 1.) Bilinear horizontal interpolation
140	! This has to be done separately for 6 fields (Temporal(2)*Vertical(3))
141	!**********************************************************************
142
143	! Loop over 2 time steps and 2 levels
144	!************************************
145
146	usl=0.
147	vsl=0.
148	wsl=0.
149	usq=0.
150	vsq=0.
151	wsq=0.
152	do m=1,2
153	indexh=memind(m)
154	do n=1,2
155	indzh=indz+n-1
156
157	u1(n)=p1*uun(ix ,jy ,indzh,indexh,ngrid) &
158	+p2*uun(ixp,jy ,indzh,indexh,ngrid) &
159	+p3*uun(ix ,jyp,indzh,indexh,ngrid) &
160	+p4*uun(ixp,jyp,indzh,indexh,ngrid)
161	v1(n)=p1*vvn(ix ,jy ,indzh,indexh,ngrid) &
162	+p2*vvn(ixp,jy ,indzh,indexh,ngrid) &
163	+p3*vvn(ix ,jyp,indzh,indexh,ngrid) &
164	+p4*vvn(ixp,jyp,indzh,indexh,ngrid)
165	w1(n)=p1*wwn(ix ,jy ,indzh,indexh,ngrid) &
166	+p2*wwn(ixp,jy ,indzh,indexh,ngrid) &
167	+p3*wwn(ix ,jyp,indzh,indexh,ngrid) &
168	+p4*wwn(ixp,jyp,indzh,indexh,ngrid)
169
170	usl=usl+uun(ix ,jy ,indzh,indexh,ngrid)+ &
171	uun(ixp,jy ,indzh,indexh,ngrid) &
172	+uun(ix ,jyp,indzh,indexh,ngrid)+ &
173	uun(ixp,jyp,indzh,indexh,ngrid)
174	vsl=vsl+vvn(ix ,jy ,indzh,indexh,ngrid)+ &
175	vvn(ixp,jy ,indzh,indexh,ngrid) &
176	+vvn(ix ,jyp,indzh,indexh,ngrid)+ &
177	vvn(ixp,jyp,indzh,indexh,ngrid)
178	wsl=wsl+wwn(ix ,jy ,indzh,indexh,ngrid)+ &
179	wwn(ixp,jy ,indzh,indexh,ngrid) &
180	+wwn(ix ,jyp,indzh,indexh,ngrid)+ &
181	wwn(ixp,jyp,indzh,indexh,ngrid)
182
183	usq=usq+uun(ix ,jy ,indzh,indexh,ngrid)* &
184	uun(ix ,jy ,indzh,indexh,ngrid)+ &
185	uun(ixp,jy ,indzh,indexh,ngrid)*uun(ixp,jy ,indzh,indexh,ngrid)+ &
186	uun(ix ,jyp,indzh,indexh,ngrid)*uun(ix ,jyp,indzh,indexh,ngrid)+ &
187	uun(ixp,jyp,indzh,indexh,ngrid)*uun(ixp,jyp,indzh,indexh,ngrid)
188	vsq=vsq+vvn(ix ,jy ,indzh,indexh,ngrid)* &
189	vvn(ix ,jy ,indzh,indexh,ngrid)+ &
190	vvn(ixp,jy ,indzh,indexh,ngrid)*vvn(ixp,jy ,indzh,indexh,ngrid)+ &
191	vvn(ix ,jyp,indzh,indexh,ngrid)*vvn(ix ,jyp,indzh,indexh,ngrid)+ &
192	vvn(ixp,jyp,indzh,indexh,ngrid)*vvn(ixp,jyp,indzh,indexh,ngrid)
193	wsq=wsq+wwn(ix ,jy ,indzh,indexh,ngrid)* &
194	wwn(ix ,jy ,indzh,indexh,ngrid)+ &
195	wwn(ixp,jy ,indzh,indexh,ngrid)*wwn(ixp,jy ,indzh,indexh,ngrid)+ &
196	wwn(ix ,jyp,indzh,indexh,ngrid)*wwn(ix ,jyp,indzh,indexh,ngrid)+ &
197	wwn(ixp,jyp,indzh,indexh,ngrid)*wwn(ixp,jyp,indzh,indexh,ngrid)
198	end do
199
200
201	!**********************************
202	! 2.) Linear vertical interpolation
203	!**********************************
204
205	uh(m)=dz2u1(1)+dz1u1(2)
206	vh(m)=dz2v1(1)+dz1v1(2)
207	wh(m)=dz2w1(1)+dz1w1(2)
208	end do
209
210
211	!************************************
212	! 3.) Temporal interpolation (linear)
213	!************************************
214
215	u=(uh(1)dt2+uh(2)dt1)*dtt
216	v=(vh(1)dt2+vh(2)dt1)*dtt
217	w=(wh(1)dt2+wh(2)dt1)*dtt
218
219
220	! Compute standard deviations
221	!****************************
222
223	xaux=usq-usl*usl/16.
224	if (xaux.lt.eps) then
225	usig=0.
226	else
227	usig=sqrt(xaux/15.)
228	endif
229
230	xaux=vsq-vsl*vsl/16.
231	if (xaux.lt.eps) then
232	vsig=0.
233	else
234	vsig=sqrt(xaux/15.)
235	endif
236
237
238	xaux=wsq-wsl*wsl/16.
239	if (xaux.lt.eps) then
240	wsig=0.
241	else
242	wsig=sqrt(xaux/15.)
243	endif
244
245	end subroutine interpol_wind_nests

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format