source: branches/jerome/src_flexwrf_v3.1/interpol_wind.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: 10.0 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_wind(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
30  !                           i   i  i  i
31  !*****************************************************************************
32  !                                                                            *
33  !  This subroutine interpolates the wind data to current trajectory position.*
34  !                                                                            *
35  !    Author: A. Stohl                                                        *
36  !                                                                            *
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
102       mu =p1*m_x(ix ,jy ,1) &
103         + p2*m_x(ixp,jy ,1) &
104         + p3*m_x(ix ,jyp,1) &
105         + p4*m_x(ixp,jyp,1)
106       mv =p1*m_y(ix ,jy ,1) &
107         + p2*m_y(ixp,jy ,1) &
108         + p3*m_y(ix ,jyp,1) &
109         + p4*m_y(ixp,jyp,1)
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
1276   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  ! 1.) Bilinear horizontal interpolation
139  ! This has to be done separately for 6 fields (Temporal(2)*Vertical(3))
140  !**********************************************************************
141
142  ! Loop over 2 time steps and 2 levels
143  !************************************
144
145  usl=0.
146  vsl=0.
147  wsl=0.
148  usq=0.
149  vsq=0.
150  wsq=0.
151  do m=1,2
152    indexh=memind(m)
153    do n=1,2
154      indzh=indz+n-1
155
156      if (ngrid.lt.0) then
157        u1(n)=p1*uupol(ix ,jy ,indzh,indexh) &
158             +p2*uupol(ixp,jy ,indzh,indexh) &
159             +p3*uupol(ix ,jyp,indzh,indexh) &
160             +p4*uupol(ixp,jyp,indzh,indexh)
161        v1(n)=p1*vvpol(ix ,jy ,indzh,indexh) &
162             +p2*vvpol(ixp,jy ,indzh,indexh) &
163             +p3*vvpol(ix ,jyp,indzh,indexh) &
164             +p4*vvpol(ixp,jyp,indzh,indexh)
165        usl=usl+uupol(ix ,jy ,indzh,indexh)+ &
166             uupol(ixp,jy ,indzh,indexh) &
167             +uupol(ix ,jyp,indzh,indexh)+uupol(ixp,jyp,indzh,indexh)
168        vsl=vsl+vvpol(ix ,jy ,indzh,indexh)+ &
169             vvpol(ixp,jy ,indzh,indexh) &
170             +vvpol(ix ,jyp,indzh,indexh)+vvpol(ixp,jyp,indzh,indexh)
171
172        usq=usq+uupol(ix ,jy ,indzh,indexh)* &
173             uupol(ix ,jy ,indzh,indexh)+ &
174             uupol(ixp,jy ,indzh,indexh)*uupol(ixp,jy ,indzh,indexh)+ &
175             uupol(ix ,jyp,indzh,indexh)*uupol(ix ,jyp,indzh,indexh)+ &
176             uupol(ixp,jyp,indzh,indexh)*uupol(ixp,jyp,indzh,indexh)
177        vsq=vsq+vvpol(ix ,jy ,indzh,indexh)* &
178             vvpol(ix ,jy ,indzh,indexh)+ &
179             vvpol(ixp,jy ,indzh,indexh)*vvpol(ixp,jy ,indzh,indexh)+ &
180             vvpol(ix ,jyp,indzh,indexh)*vvpol(ix ,jyp,indzh,indexh)+ &
181             vvpol(ixp,jyp,indzh,indexh)*vvpol(ixp,jyp,indzh,indexh)
182      else
183        u1(n)=p1*uu(ix ,jy ,indzh,indexh) &
184             +p2*uu(ixp,jy ,indzh,indexh) &
185             +p3*uu(ix ,jyp,indzh,indexh) &
186             +p4*uu(ixp,jyp,indzh,indexh)
187        v1(n)=p1*vv(ix ,jy ,indzh,indexh) &
188             +p2*vv(ixp,jy ,indzh,indexh) &
189             +p3*vv(ix ,jyp,indzh,indexh) &
190             +p4*vv(ixp,jyp,indzh,indexh)
191        usl=usl+uu(ix ,jy ,indzh,indexh)+uu(ixp,jy ,indzh,indexh) &
192             +uu(ix ,jyp,indzh,indexh)+uu(ixp,jyp,indzh,indexh)
193        vsl=vsl+vv(ix ,jy ,indzh,indexh)+vv(ixp,jy ,indzh,indexh) &
194             +vv(ix ,jyp,indzh,indexh)+vv(ixp,jyp,indzh,indexh)
195
196        usq=usq+uu(ix ,jy ,indzh,indexh)*uu(ix ,jy ,indzh,indexh)+ &
197             uu(ixp,jy ,indzh,indexh)*uu(ixp,jy ,indzh,indexh)+ &
198             uu(ix ,jyp,indzh,indexh)*uu(ix ,jyp,indzh,indexh)+ &
199             uu(ixp,jyp,indzh,indexh)*uu(ixp,jyp,indzh,indexh)
200        vsq=vsq+vv(ix ,jy ,indzh,indexh)*vv(ix ,jy ,indzh,indexh)+ &
201             vv(ixp,jy ,indzh,indexh)*vv(ixp,jy ,indzh,indexh)+ &
202             vv(ix ,jyp,indzh,indexh)*vv(ix ,jyp,indzh,indexh)+ &
203             vv(ixp,jyp,indzh,indexh)*vv(ixp,jyp,indzh,indexh)
204      endif
205      w1(n)=p1*ww(ix ,jy ,indzh,indexh) &
206           +p2*ww(ixp,jy ,indzh,indexh) &
207           +p3*ww(ix ,jyp,indzh,indexh) &
208           +p4*ww(ixp,jyp,indzh,indexh)
209      wsl=wsl+ww(ix ,jy ,indzh,indexh)+ww(ixp,jy ,indzh,indexh) &
210           +ww(ix ,jyp,indzh,indexh)+ww(ixp,jyp,indzh,indexh)
211      wsq=wsq+ww(ix ,jy ,indzh,indexh)*ww(ix ,jy ,indzh,indexh)+ &
212           ww(ixp,jy ,indzh,indexh)*ww(ixp,jy ,indzh,indexh)+ &
213           ww(ix ,jyp,indzh,indexh)*ww(ix ,jyp,indzh,indexh)+ &
214           ww(ixp,jyp,indzh,indexh)*ww(ixp,jyp,indzh,indexh)
215    end do
216
217
218  !**********************************
219  ! 2.) Linear vertical interpolation
220  !**********************************
221
222    uh(m)=dz2*u1(1)+dz1*u1(2)
223    vh(m)=dz2*v1(1)+dz1*v1(2)
224    wh(m)=dz2*w1(1)+dz1*w1(2)
225  end do
226
227
228  !************************************
229  ! 3.) Temporal interpolation (linear)
230  !************************************
231
232  u=(uh(1)*dt2+uh(2)*dt1)*dtt
233  v=(vh(1)*dt2+vh(2)*dt1)*dtt
234  w=(wh(1)*dt2+wh(2)*dt1)*dtt
235
236
237  ! Compute standard deviations
238  !****************************
239
240  xaux=usq-usl*usl/16.
241  if (xaux.lt.eps) then
242    usig=0.
243  else
244    usig=sqrt(xaux/15.)
245  endif
246  if (usig.gt.1.e6) print*,'pb xaux',usq,usl
247
248  xaux=vsq-vsl*vsl/16.
249  if (xaux.lt.eps) then
250    vsig=0.
251  else
252    vsig=sqrt(xaux/15.)
253  endif
254
255
256  xaux=wsq-wsl*wsl/16.
257  if (xaux.lt.eps) then
258    wsig=0.
259  else
260    wsig=sqrt(xaux/15.)
261  endif
262
263end subroutine interpol_wind
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG