Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

interpol_all_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: 11.3 KB

Line
1	!***********************************************************************
2	!* Copyright 2012,2013 *
3	!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine, *
4	!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5	!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso, *
6	!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann, *
7	!* *
8	!* This file is part of FLEXPART WRF *
9	!* *
10	!* FLEXPART is free software: you can redistribute it and/or modify *
11	!* it under the terms of the GNU General Public License as published by*
12	!* the Free Software Foundation, either version 3 of the License, or *
13	!* (at your option) any later version. *
14	!* *
15	!* FLEXPART is distributed in the hope that it will be useful, *
16	!* but WITHOUT ANY WARRANTY; without even the implied warranty of *
17	!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18	!* GNU General Public License for more details. *
19	!* *
20	!* You should have received a copy of the GNU General Public License *
21	!* along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
22	!***********************************************************************
23	subroutine interpol_all_nests(itime,xt,yt,zt, &
24	uprof,vprof,wprof, usigprof,vsigprof,wsigprof, &
25	rhoprof,rhogradprof, tkeprof,pttprof, &
26	u,v,w,usig,vsig,wsig,pvi, &
27	p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2, &
28	ix,jy,ixp,jyp,ngrid,indz,indzp, depoindicator, &
29	indzindicator, &
30	ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
31	sigw,dsigwdz,dsigw2dz,mu,mv)
32
33	! i i i i
34	!*******************************************************************************
35	! *
36	! This subroutine interpolates everything that is needed for calculating the *
37	! dispersion. *
38	! Version for interpolating nested grids. *
39	! *
40	! Author: A. Stohl *
41	! *
42	! 9 February 1999 *
43	! 16 December 1997 *
44	! *
45	! Revision March 2005 by AST : all output variables in common block *
46	! calculation of standard deviation done in this *
47	! routine rather than subroutine call in order *
48	! to save computation time *
49	! *
50	!*******************************************************************************
51	! *
52	! Variables: *
53	! itime [s] current temporal position *
54	! memtime(3) [s] times of the wind fields in memory *
55	! xt,yt,zt coordinates position for which wind data shall be calculat*
56	! *
57	! Constants: *
58	! *
59	!*******************************************************************************
60	! 12 JUNE 2007 compute tkeprof, add a variable y4(2)
61	! 25 June 2007 compute pttprof, y5
62	! compute tkeprof for all levels
63	!*******************************************************************************
64	use par_mod
65	use com_mod
66	! use interpol_mod
67	! use hanna_mod
68
69	implicit none
70
71	integer :: itime
72	real :: xt,yt,zt
73
74	! Auxiliary variables needed for interpolation
75	real :: ust1(2),wst1(2),oli1(2),oliaux
76	real :: y1(2),y2(2),y3(2),rho1(2),rhograd1(2),y4(2),y5(2)
77	real :: usl,vsl,wsl,usq,vsq,wsq,xaux
78	integer :: i,m,n,indexh
79	real,parameter :: eps=1.0e-30
80	real :: uprof(nzmax),vprof(nzmax),wprof(nzmax)
81	real :: usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax)
82	real :: rhoprof(nzmax),rhogradprof(nzmax)
83	real :: tkeprof(nzmax),pttprof(nzmax)
84	real :: u,v,w,usig,vsig,wsig,pvi,mu,mv
85
86	real :: p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2
87	integer :: ix,jy,ixp,jyp,ngrid,indz,indzp
88	logical :: depoindicator(maxspec)
89	logical :: indzindicator(nzmax)
90	real :: ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw
91	real :: sigw,dsigwdz,dsigw2dz
92
93
94	!********************************************
95	! Multilinear interpolation in time and space
96	!********************************************
97
98	! Determine the lower left corner and its distance to the current position
99	!*************************************************************************
100
101	ddx=xt-real(ix)
102	ddy=yt-real(jy)
103	rddx=1.-ddx
104	rddy=1.-ddy
105	p1=rddx*rddy
106	p2=ddx*rddy
107	p3=rddx*ddy
108	p4=ddx*ddy
109
110	! Calculate variables for time interpolation
111	!*******************************************
112
113	dt1=real(itime-memtime(1))
114	dt2=real(memtime(2)-itime)
115	dtt=1./(dt1+dt2)
116
117
118	!*****************************************
119	! 1. Interpolate u, w and Obukhov length
120	!*****************************************
121
122	! a) Bilinear horizontal interpolation
123	do m=1,2
124	indexh=memind(m)
125
126	ust1(m)=p1*ustarn(ix ,jy ,1,indexh,ngrid) &
127	+ p2*ustarn(ixp,jy ,1,indexh,ngrid) &
128	+ p3*ustarn(ix ,jyp,1,indexh,ngrid) &
129	+ p4*ustarn(ixp,jyp,1,indexh,ngrid)
130	wst1(m)=p1*wstarn(ix ,jy ,1,indexh,ngrid) &
131	+ p2*wstarn(ixp,jy ,1,indexh,ngrid) &
132	+ p3*wstarn(ix ,jyp,1,indexh,ngrid) &
133	+ p4*wstarn(ixp,jyp,1,indexh,ngrid)
134	oli1(m)=p1*olin(ix ,jy ,1,indexh,ngrid) &
135	+ p2*olin(ixp,jy ,1,indexh,ngrid) &
136	+ p3*olin(ix ,jyp,1,indexh,ngrid) &
137	+ p4*olin(ixp,jyp,1,indexh,ngrid)
138	end do
139	mu =p1*m_xn(ix ,jy ,1,ngrid) &
140	+ p2*m_xn(ixp,jy ,1,ngrid) &
141	+ p3*m_xn(ix ,jyp,1,ngrid) &
142	+ p4*m_xn(ixp,jyp,1,ngrid)
143	mv =p1*m_yn(ix ,jy ,1,ngrid) &
144	+ p2*m_yn(ixp,jy ,1,ngrid) &
145	+ p3*m_yn(ix ,jyp,1,ngrid) &
146	+ p4*m_yn(ixp,jyp,1,ngrid)
147
148
149	! b) Temporal interpolation
150
151	ust=(ust1(1)dt2+ust1(2)dt1)*dtt
152	wst=(wst1(1)dt2+wst1(2)dt1)*dtt
153	oliaux=(oli1(1)dt2+oli1(2)dt1)*dtt
154
155	if (oliaux.ne.0.) then
156	ol=1./oliaux
157	else
158	ol=99999.
159	endif
160
161
162	!*****************************************************
163	! 2. Interpolate vertical profiles of u,v,w,rho,drhodz
164	!*****************************************************
165
166
167	! Determine the level below the current position
168	!***********************************************
169
170	do i=2,nz
171	if (height(i).gt.zt) then
172	indz=i-1
173	indzp=i
174	goto 6
175	endif
176	end do
177	6 continue
178
179	!**************************************
180	! 1.) Bilinear horizontal interpolation
181	! 2.) Temporal interpolation (linear)
182	!**************************************
183
184	! Loop over 2 time steps and indz levels
185	!***************************************
186
187	do n=indz,indz+1
188	usl=0.
189	vsl=0.
190	wsl=0.
191	usq=0.
192	vsq=0.
193	wsq=0.
194	do m=1,2
195	indexh=memind(m)
196	y1(m)=p1*uun(ix ,jy ,n,indexh,ngrid) &
197	+p2*uun(ixp,jy ,n,indexh,ngrid) &
198	+p3*uun(ix ,jyp,n,indexh,ngrid) &
199	+p4*uun(ixp,jyp,n,indexh,ngrid)
200	y2(m)=p1*vvn(ix ,jy ,n,indexh,ngrid) &
201	+p2*vvn(ixp,jy ,n,indexh,ngrid) &
202	+p3*vvn(ix ,jyp,n,indexh,ngrid) &
203	+p4*vvn(ixp,jyp,n,indexh,ngrid)
204	y3(m)=p1*wwn(ix ,jy ,n,indexh,ngrid) &
205	+p2*wwn(ixp,jy ,n,indexh,ngrid) &
206	+p3*wwn(ix ,jyp,n,indexh,ngrid) &
207	+p4*wwn(ixp,jyp,n,indexh,ngrid)
208	rhograd1(m)=p1*drhodzn(ix ,jy ,n,indexh,ngrid) &
209	+p2*drhodzn(ixp,jy ,n,indexh,ngrid) &
210	+p3*drhodzn(ix ,jyp,n,indexh,ngrid) &
211	+p4*drhodzn(ixp,jyp,n,indexh,ngrid)
212	rho1(m)=p1*rhon(ix ,jy ,n,indexh,ngrid) &
213	+p2*rhon(ixp,jy ,n,indexh,ngrid) &
214	+p3*rhon(ix ,jyp,n,indexh,ngrid) &
215	+p4*rhon(ixp,jyp,n,indexh,ngrid)
216
217
218	usl=usl+uun(ix ,jy ,n,indexh,ngrid)+uun(ixp,jy ,n,indexh,ngrid) &
219	+uun(ix ,jyp,n,indexh,ngrid)+uun(ixp,jyp,n,indexh,ngrid)
220	vsl=vsl+vvn(ix ,jy ,n,indexh,ngrid)+vvn(ixp,jy ,n,indexh,ngrid) &
221	+vvn(ix ,jyp,n,indexh,ngrid)+vvn(ixp,jyp,n,indexh,ngrid)
222	wsl=wsl+wwn(ix ,jy ,n,indexh,ngrid)+wwn(ixp,jy ,n,indexh,ngrid) &
223	+wwn(ix ,jyp,n,indexh,ngrid)+wwn(ixp,jyp,n,indexh,ngrid)
224
225	usq=usq+uun(ix ,jy ,n,indexh,ngrid)*uun(ix ,jy ,n,indexh,ngrid)+ &
226	uun(ixp,jy ,n,indexh,ngrid)*uun(ixp,jy ,n,indexh,ngrid)+ &
227	uun(ix ,jyp,n,indexh,ngrid)*uun(ix ,jyp,n,indexh,ngrid)+ &
228	uun(ixp,jyp,n,indexh,ngrid)*uun(ixp,jyp,n,indexh,ngrid)
229	vsq=vsq+vvn(ix ,jy ,n,indexh,ngrid)*vvn(ix ,jy ,n,indexh,ngrid)+ &
230	vvn(ixp,jy ,n,indexh,ngrid)*vvn(ixp,jy ,n,indexh,ngrid)+ &
231	vvn(ix ,jyp,n,indexh,ngrid)*vvn(ix ,jyp,n,indexh,ngrid)+ &
232	vvn(ixp,jyp,n,indexh,ngrid)*vvn(ixp,jyp,n,indexh,ngrid)
233	wsq=wsq+wwn(ix ,jy ,n,indexh,ngrid)*wwn(ix ,jy ,n,indexh,ngrid)+ &
234	wwn(ixp,jy ,n,indexh,ngrid)*wwn(ixp,jy ,n,indexh,ngrid)+ &
235	wwn(ix ,jyp,n,indexh,ngrid)*wwn(ix ,jyp,n,indexh,ngrid)+ &
236	wwn(ixp,jyp,n,indexh,ngrid)*wwn(ixp,jyp,n,indexh,ngrid)
237	end do
238	uprof(n)=(y1(1)dt2+y1(2)dt1)*dtt
239	vprof(n)=(y2(1)dt2+y2(2)dt1)*dtt
240	wprof(n)=(y3(1)dt2+y3(2)dt1)*dtt
241	rhoprof(n)=(rho1(1)dt2+rho1(2)dt1)*dtt
242	rhogradprof(n)=(rhograd1(1)dt2+rhograd1(2)dt1)*dtt
243	indzindicator(n)=.false.
244
245	! Compute standard deviations
246	!****************************
247
248	xaux=usq-usl*usl/8.
249	if (xaux.lt.eps) then
250	usigprof(n)=0.
251	else
252	usigprof(n)=sqrt(xaux/7.)
253	endif
254
255	xaux=vsq-vsl*vsl/8.
256	if (xaux.lt.eps) then
257	vsigprof(n)=0.
258	else
259	vsigprof(n)=sqrt(xaux/7.)
260	endif
261
262
263	xaux=wsq-wsl*wsl/8.
264	if (xaux.lt.eps) then
265	wsigprof(n)=0.
266	else
267	wsigprof(n)=sqrt(xaux/7.)
268	endif
269
270	end do
271
272	! compute tke profile for all levels
273
274	do n=1,nz
275	do m=1,2
276	indexh=memind(m)
277	y4(m)=p1*tken(ix ,jy ,n,indexh,ngrid) &
278	+p2*tken(ixp,jy ,n,indexh,ngrid) &
279	+p3*tken(ix ,jyp,n,indexh,ngrid) &
280	+p4*tken(ixp,jyp,n,indexh,ngrid)
281	y5(m)=p1*pttn(ix ,jy ,n,indexh,ngrid) &
282	+p2*pttn(ixp,jy ,n,indexh,ngrid) &
283	+p3*pttn(ix ,jyp,n,indexh,ngrid) &
284	+p4*pttn(ixp,jyp,n,indexh,ngrid)
285
286	enddo
287	tkeprof(n)=(y4(1)dt2+y4(2)dt1)*dtt
288	pttprof(n)=(y5(1)dt2+y5(2)dt1)*dtt
289	enddo
290
291
292	end subroutine interpol_all_nests
293

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

Download in other formats:

Original Format