CTBTO FLEXPART WO4 (2015-10-15)
 All Classes Files Functions Variables
interpol_all.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_all(itime,xt,yt,zt)
23  ! i i i i
24  !*****************************************************************************
25  ! *
26  ! This subroutine interpolates everything that is needed for calculating the*
27  ! dispersion. *
28  ! *
29  ! Author: A. Stohl *
30  ! *
31  ! 16 December 1997 *
32  ! *
33  ! Revision March 2005 by AST : all output variables in common block cal- *
34  ! culation of standard deviation done in this *
35  ! routine rather than subroutine call in order *
36  ! to save computation time *
37  ! *
38  !*****************************************************************************
39  ! *
40  ! Variables: *
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  ! culated *
45  ! *
46  ! Constants: *
47  ! *
48  !*****************************************************************************
49 
50  use par_mod
51  use com_mod
52  use interpol_mod
53  use hanna_mod
54 
55  implicit none
56 
57  integer :: itime
58  real :: xt,yt,zt
59 
60  ! Auxiliary variables needed for interpolation
61  real :: ust1(2),wst1(2),oli1(2),oliaux
62  real :: y1(2),y2(2),y3(2),rho1(2),rhograd1(2)
63  real :: usl,vsl,wsl,usq,vsq,wsq,xaux
64  integer :: i,m,n,indexh
65  real,parameter :: eps=1.0e-30
66 
67 
68  !********************************************
69  ! Multilinear interpolation in time and space
70  !********************************************
71 
72  ! Determine the lower left corner and its distance to the current position
73  !*************************************************************************
74 
75  ddx=xt-real(ix)
76  ddy=yt-real(jy)
77  rddx=1.-ddx
78  rddy=1.-ddy
79  p1=rddx*rddy
80  p2=ddx*rddy
81  p3=rddx*ddy
82  p4=ddx*ddy
83 
84  ! Calculate variables for time interpolation
85  !*******************************************
86 
87  dt1=real(itime-memtime(1))
88  dt2=real(memtime(2)-itime)
89  dtt=1./(dt1+dt2)
90 
91 
92  !*****************************************
93  ! 1. Interpolate u*, w* and Obukhov length
94  !*****************************************
95 
96  ! a) Bilinear horizontal interpolation
97 
98  do m=1,2
99  indexh=memind(m)
100 
101  ust1(m)=p1*ustar(ix ,jy ,1,indexh) &
102  + p2*ustar(ixp,jy ,1,indexh) &
103  + p3*ustar(ix ,jyp,1,indexh) &
104  + p4*ustar(ixp,jyp,1,indexh)
105  wst1(m)=p1*wstar(ix ,jy ,1,indexh) &
106  + p2*wstar(ixp,jy ,1,indexh) &
107  + p3*wstar(ix ,jyp,1,indexh) &
108  + p4*wstar(ixp,jyp,1,indexh)
109  oli1(m)=p1*oli(ix ,jy ,1,indexh) &
110  + p2*oli(ixp,jy ,1,indexh) &
111  + p3*oli(ix ,jyp,1,indexh) &
112  + p4*oli(ixp,jyp,1,indexh)
113  end do
114 
115  ! b) Temporal interpolation
116 
117  ust=(ust1(1)*dt2+ust1(2)*dt1)*dtt
118  wst=(wst1(1)*dt2+wst1(2)*dt1)*dtt
119  oliaux=(oli1(1)*dt2+oli1(2)*dt1)*dtt
120 
121  if (oliaux.ne.0.) then
122  ol=1./oliaux
123  else
124  ol=99999.
125  endif
126 
127 
128  !*****************************************************
129  ! 2. Interpolate vertical profiles of u,v,w,rho,drhodz
130  !*****************************************************
131 
132 
133  ! Determine the level below the current position
134  !***********************************************
135 
136  do i=2,nz
137  if (height(i).gt.zt) then
138  indz=i-1
139  indzp=i
140  goto 6
141  endif
142  end do
143 6 continue
144 
145  !**************************************
146  ! 1.) Bilinear horizontal interpolation
147  ! 2.) Temporal interpolation (linear)
148  !**************************************
149 
150  ! Loop over 2 time steps and indz levels
151  !***************************************
152 
153  do n=indz,indzp
154  usl=0.
155  vsl=0.
156  wsl=0.
157  usq=0.
158  vsq=0.
159  wsq=0.
160  do m=1,2
161  indexh=memind(m)
162  if (ngrid.lt.0) then
163  y1(m)=p1*uupol(ix ,jy ,n,indexh) &
164  +p2*uupol(ixp,jy ,n,indexh) &
165  +p3*uupol(ix ,jyp,n,indexh) &
166  +p4*uupol(ixp,jyp,n,indexh)
167  y2(m)=p1*vvpol(ix ,jy ,n,indexh) &
168  +p2*vvpol(ixp,jy ,n,indexh) &
169  +p3*vvpol(ix ,jyp,n,indexh) &
170  +p4*vvpol(ixp,jyp,n,indexh)
171  usl=usl+uupol(ix ,jy ,n,indexh)+uupol(ixp,jy ,n,indexh) &
172  +uupol(ix ,jyp,n,indexh)+uupol(ixp,jyp,n,indexh)
173  vsl=vsl+vvpol(ix ,jy ,n,indexh)+vvpol(ixp,jy ,n,indexh) &
174  +vvpol(ix ,jyp,n,indexh)+vvpol(ixp,jyp,n,indexh)
175 
176  usq=usq+uupol(ix ,jy ,n,indexh)*uupol(ix ,jy ,n,indexh)+ &
177  uupol(ixp,jy ,n,indexh)*uupol(ixp,jy ,n,indexh)+ &
178  uupol(ix ,jyp,n,indexh)*uupol(ix ,jyp,n,indexh)+ &
179  uupol(ixp,jyp,n,indexh)*uupol(ixp,jyp,n,indexh)
180  vsq=vsq+vvpol(ix ,jy ,n,indexh)*vvpol(ix ,jy ,n,indexh)+ &
181  vvpol(ixp,jy ,n,indexh)*vvpol(ixp,jy ,n,indexh)+ &
182  vvpol(ix ,jyp,n,indexh)*vvpol(ix ,jyp,n,indexh)+ &
183  vvpol(ixp,jyp,n,indexh)*vvpol(ixp,jyp,n,indexh)
184  else
185  y1(m)=p1*uu(ix ,jy ,n,indexh) &
186  +p2*uu(ixp,jy ,n,indexh) &
187  +p3*uu(ix ,jyp,n,indexh) &
188  +p4*uu(ixp,jyp,n,indexh)
189  y2(m)=p1*vv(ix ,jy ,n,indexh) &
190  +p2*vv(ixp,jy ,n,indexh) &
191  +p3*vv(ix ,jyp,n,indexh) &
192  +p4*vv(ixp,jyp,n,indexh)
193  usl=usl+uu(ix ,jy ,n,indexh)+uu(ixp,jy ,n,indexh) &
194  +uu(ix ,jyp,n,indexh)+uu(ixp,jyp,n,indexh)
195  vsl=vsl+vv(ix ,jy ,n,indexh)+vv(ixp,jy ,n,indexh) &
196  +vv(ix ,jyp,n,indexh)+vv(ixp,jyp,n,indexh)
197 
198  usq=usq+uu(ix ,jy ,n,indexh)*uu(ix ,jy ,n,indexh)+ &
199  uu(ixp,jy ,n,indexh)*uu(ixp,jy ,n,indexh)+ &
200  uu(ix ,jyp,n,indexh)*uu(ix ,jyp,n,indexh)+ &
201  uu(ixp,jyp,n,indexh)*uu(ixp,jyp,n,indexh)
202  vsq=vsq+vv(ix ,jy ,n,indexh)*vv(ix ,jy ,n,indexh)+ &
203  vv(ixp,jy ,n,indexh)*vv(ixp,jy ,n,indexh)+ &
204  vv(ix ,jyp,n,indexh)*vv(ix ,jyp,n,indexh)+ &
205  vv(ixp,jyp,n,indexh)*vv(ixp,jyp,n,indexh)
206  endif
207  y3(m)=p1*ww(ix ,jy ,n,indexh) &
208  +p2*ww(ixp,jy ,n,indexh) &
209  +p3*ww(ix ,jyp,n,indexh) &
210  +p4*ww(ixp,jyp,n,indexh)
211  rhograd1(m)=p1*drhodz(ix ,jy ,n,indexh) &
212  +p2*drhodz(ixp,jy ,n,indexh) &
213  +p3*drhodz(ix ,jyp,n,indexh) &
214  +p4*drhodz(ixp,jyp,n,indexh)
215  rho1(m)=p1*rho(ix ,jy ,n,indexh) &
216  +p2*rho(ixp,jy ,n,indexh) &
217  +p3*rho(ix ,jyp,n,indexh) &
218  +p4*rho(ixp,jyp,n,indexh)
219  wsl=wsl+ww(ix ,jy ,n,indexh)+ww(ixp,jy ,n,indexh) &
220  +ww(ix ,jyp,n,indexh)+ww(ixp,jyp,n,indexh)
221  wsq=wsq+ww(ix ,jy ,n,indexh)*ww(ix ,jy ,n,indexh)+ &
222  ww(ixp,jy ,n,indexh)*ww(ixp,jy ,n,indexh)+ &
223  ww(ix ,jyp,n,indexh)*ww(ix ,jyp,n,indexh)+ &
224  ww(ixp,jyp,n,indexh)*ww(ixp,jyp,n,indexh)
225  end do
226  uprof(n)=(y1(1)*dt2+y1(2)*dt1)*dtt
227  vprof(n)=(y2(1)*dt2+y2(2)*dt1)*dtt
228  wprof(n)=(y3(1)*dt2+y3(2)*dt1)*dtt
229  rhoprof(n)=(rho1(1)*dt2+rho1(2)*dt1)*dtt
230  rhogradprof(n)=(rhograd1(1)*dt2+rhograd1(2)*dt1)*dtt
231  indzindicator(n)=.false.
232 
233  ! Compute standard deviations
234  !****************************
235 
236  xaux=usq-usl*usl/8.
237  if (xaux.lt.eps) then
238  usigprof(n)=0.
239  else
240  usigprof(n)=sqrt(xaux/7.)
241  endif
242 
243  xaux=vsq-vsl*vsl/8.
244  if (xaux.lt.eps) then
245  vsigprof(n)=0.
246  else
247  vsigprof(n)=sqrt(xaux/7.)
248  endif
249 
250 
251  xaux=wsq-wsl*wsl/8.
252  if (xaux.lt.eps) then
253  wsigprof(n)=0.
254  else
255  wsigprof(n)=sqrt(xaux/7.)
256  endif
257 
258  end do
259 
260 
261 end subroutine interpol_all
interpol_all
subroutine interpol_all(itime, xt, yt, zt)
Definition: interpol_all.f90:22
par_mod
Definition: par_mod.F90:34
interpol_mod
Definition: interpol_mod.f90:1
hanna_mod
Definition: hanna_mod.f90:1
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