source: trunk/src/interpol_all.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 11 years ago
File size: 9.7 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
22subroutine 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
1436   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
261end subroutine interpol_all
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG