CTBTO FLEXPART WO4 (2015-10-15)
 All Classes Files Functions Variables
interpol_wind.f90
Go to the documentation of this file.
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(itime,xt,yt,zt)
23  ! i i i i
24  !*****************************************************************************
25  ! *
26  ! This subroutine interpolates the wind data to current trajectory position.*
27  ! *
28  ! Author: A. Stohl *
29  ! *
30  ! 16 December 1997 *
31  ! *
32  ! Revision March 2005 by AST : all output variables in common block cal- *
33  ! culation of standard deviation done in this *
34  ! routine rather than subroutine call in order *
35  ! to save computation time *
36  ! *
37  !*****************************************************************************
38  ! *
39  ! Variables: *
40  ! u,v,w wind components *
41  ! itime [s] current temporal position *
42  ! memtime(3) [s] times of the wind fields in memory *
43  ! xt,yt,zt coordinates position for which wind data shall be *
44  ! calculated *
45  ! *
46  ! Constants: *
47  ! *
48  !*****************************************************************************
49 
50  use par_mod
51  use com_mod
52  use interpol_mod
53 
54  implicit none
55 
56  integer :: itime
57  real :: xt,yt,zt
58 
59  ! Auxiliary variables needed for interpolation
60  real :: dz1,dz2,dz
61  real :: u1(2),v1(2),w1(2),uh(2),vh(2),wh(2)
62  real :: usl,vsl,wsl,usq,vsq,wsq,xaux
63  integer :: i,m,n,indexh,indzh
64  real,parameter :: eps=1.0e-30
65 
66 
67  !********************************************
68  ! Multilinear interpolation in time and space
69  !********************************************
70 
71  ! Determine the lower left corner and its distance to the current position
72  !*************************************************************************
73 
74  ddx=xt-real(ix)
75  ddy=yt-real(jy)
76  rddx=1.-ddx
77  rddy=1.-ddy
78  p1=rddx*rddy
79  p2=ddx*rddy
80  p3=rddx*ddy
81  p4=ddx*ddy
82 
83  ! Calculate variables for time interpolation
84  !*******************************************
85 
86  dt1=real(itime-memtime(1))
87  dt2=real(memtime(2)-itime)
88  dtt=1./(dt1+dt2)
89 
90  ! Determine the level below the current position for u,v
91  !*******************************************************
92 
93  do i=2,nz
94  if (height(i).gt.zt) then
95  indz=i-1
96  goto 6
97  endif
98  end do
99 6 continue
100 
101 
102  ! Vertical distance to the level below and above current position
103  !****************************************************************
104 
105  dz=1./(height(indz+1)-height(indz))
106  dz1=(zt-height(indz))*dz
107  dz2=(height(indz+1)-zt)*dz
108 
109  !**********************************************************************
110  ! 1.) Bilinear horizontal interpolation
111  ! This has to be done separately for 6 fields (Temporal(2)*Vertical(3))
112  !**********************************************************************
113 
114  ! Loop over 2 time steps and 2 levels
115  !************************************
116 
117  usl=0.
118  vsl=0.
119  wsl=0.
120  usq=0.
121  vsq=0.
122  wsq=0.
123  do m=1,2
124  indexh=memind(m)
125  do n=1,2
126  indzh=indz+n-1
127 
128  if (ngrid.lt.0) then
129  u1(n)=p1*uupol(ix ,jy ,indzh,indexh) &
130  +p2*uupol(ixp,jy ,indzh,indexh) &
131  +p3*uupol(ix ,jyp,indzh,indexh) &
132  +p4*uupol(ixp,jyp,indzh,indexh)
133  v1(n)=p1*vvpol(ix ,jy ,indzh,indexh) &
134  +p2*vvpol(ixp,jy ,indzh,indexh) &
135  +p3*vvpol(ix ,jyp,indzh,indexh) &
136  +p4*vvpol(ixp,jyp,indzh,indexh)
137  usl=usl+uupol(ix ,jy ,indzh,indexh)+ &
138  uupol(ixp,jy ,indzh,indexh) &
139  +uupol(ix ,jyp,indzh,indexh)+uupol(ixp,jyp,indzh,indexh)
140  vsl=vsl+vvpol(ix ,jy ,indzh,indexh)+ &
141  vvpol(ixp,jy ,indzh,indexh) &
142  +vvpol(ix ,jyp,indzh,indexh)+vvpol(ixp,jyp,indzh,indexh)
143 
144  usq=usq+uupol(ix ,jy ,indzh,indexh)* &
145  uupol(ix ,jy ,indzh,indexh)+ &
146  uupol(ixp,jy ,indzh,indexh)*uupol(ixp,jy ,indzh,indexh)+ &
147  uupol(ix ,jyp,indzh,indexh)*uupol(ix ,jyp,indzh,indexh)+ &
148  uupol(ixp,jyp,indzh,indexh)*uupol(ixp,jyp,indzh,indexh)
149  vsq=vsq+vvpol(ix ,jy ,indzh,indexh)* &
150  vvpol(ix ,jy ,indzh,indexh)+ &
151  vvpol(ixp,jy ,indzh,indexh)*vvpol(ixp,jy ,indzh,indexh)+ &
152  vvpol(ix ,jyp,indzh,indexh)*vvpol(ix ,jyp,indzh,indexh)+ &
153  vvpol(ixp,jyp,indzh,indexh)*vvpol(ixp,jyp,indzh,indexh)
154  else
155  u1(n)=p1*uu(ix ,jy ,indzh,indexh) &
156  +p2*uu(ixp,jy ,indzh,indexh) &
157  +p3*uu(ix ,jyp,indzh,indexh) &
158  +p4*uu(ixp,jyp,indzh,indexh)
159  v1(n)=p1*vv(ix ,jy ,indzh,indexh) &
160  +p2*vv(ixp,jy ,indzh,indexh) &
161  +p3*vv(ix ,jyp,indzh,indexh) &
162  +p4*vv(ixp,jyp,indzh,indexh)
163  usl=usl+uu(ix ,jy ,indzh,indexh)+uu(ixp,jy ,indzh,indexh) &
164  +uu(ix ,jyp,indzh,indexh)+uu(ixp,jyp,indzh,indexh)
165  vsl=vsl+vv(ix ,jy ,indzh,indexh)+vv(ixp,jy ,indzh,indexh) &
166  +vv(ix ,jyp,indzh,indexh)+vv(ixp,jyp,indzh,indexh)
167 
168  usq=usq+uu(ix ,jy ,indzh,indexh)*uu(ix ,jy ,indzh,indexh)+ &
169  uu(ixp,jy ,indzh,indexh)*uu(ixp,jy ,indzh,indexh)+ &
170  uu(ix ,jyp,indzh,indexh)*uu(ix ,jyp,indzh,indexh)+ &
171  uu(ixp,jyp,indzh,indexh)*uu(ixp,jyp,indzh,indexh)
172  vsq=vsq+vv(ix ,jy ,indzh,indexh)*vv(ix ,jy ,indzh,indexh)+ &
173  vv(ixp,jy ,indzh,indexh)*vv(ixp,jy ,indzh,indexh)+ &
174  vv(ix ,jyp,indzh,indexh)*vv(ix ,jyp,indzh,indexh)+ &
175  vv(ixp,jyp,indzh,indexh)*vv(ixp,jyp,indzh,indexh)
176  endif
177  w1(n)=p1*ww(ix ,jy ,indzh,indexh) &
178  +p2*ww(ixp,jy ,indzh,indexh) &
179  +p3*ww(ix ,jyp,indzh,indexh) &
180  +p4*ww(ixp,jyp,indzh,indexh)
181  wsl=wsl+ww(ix ,jy ,indzh,indexh)+ww(ixp,jy ,indzh,indexh) &
182  +ww(ix ,jyp,indzh,indexh)+ww(ixp,jyp,indzh,indexh)
183  wsq=wsq+ww(ix ,jy ,indzh,indexh)*ww(ix ,jy ,indzh,indexh)+ &
184  ww(ixp,jy ,indzh,indexh)*ww(ixp,jy ,indzh,indexh)+ &
185  ww(ix ,jyp,indzh,indexh)*ww(ix ,jyp,indzh,indexh)+ &
186  ww(ixp,jyp,indzh,indexh)*ww(ixp,jyp,indzh,indexh)
187  end do
188 
189 
190  !**********************************
191  ! 2.) Linear vertical interpolation
192  !**********************************
193 
194  uh(m)=dz2*u1(1)+dz1*u1(2)
195  vh(m)=dz2*v1(1)+dz1*v1(2)
196  wh(m)=dz2*w1(1)+dz1*w1(2)
197  end do
198 
199 
200  !************************************
201  ! 3.) Temporal interpolation (linear)
202  !************************************
203 
204  u=(uh(1)*dt2+uh(2)*dt1)*dtt
205  v=(vh(1)*dt2+vh(2)*dt1)*dtt
206  w=(wh(1)*dt2+wh(2)*dt1)*dtt
207 
208 
209  ! Compute standard deviations
210  !****************************
211 
212  xaux=usq-usl*usl/16.
213  if (xaux.lt.eps) then
214  usig=0.
215  else
216  usig=sqrt(xaux/15.)
217  endif
218 
219  xaux=vsq-vsl*vsl/16.
220  if (xaux.lt.eps) then
221  vsig=0.
222  else
223  vsig=sqrt(xaux/15.)
224  endif
225 
226 
227  xaux=wsq-wsl*wsl/16.
228  if (xaux.lt.eps) then
229  wsig=0.
230  else
231  wsig=sqrt(xaux/15.)
232  endif
233 
234 end subroutine interpol_wind
par_mod
Definition: par_mod.F90:34
interpol_mod
Definition: interpol_mod.f90:1
interpol_wind
subroutine interpol_wind(itime, xt, yt, zt)
Definition: interpol_wind.f90:22
com_mod
Definition: com_mod.f90:13
Generated on Thu Oct 15 2015 06:20:58 for CTBTO FLEXPART WO4 (2015-10-15) by   doxygen 1.8.6