source: flexpart.git/src/interpol_wind.f90 @ 02095e3

FPv9.3.1FPv9.3.1b_testingFPv9.3.2NetCDFdepositiondevflexpart-noresmfp9.3.1-20161214-nc4grib2nc4_repairinputlistlaptoprelease-10svn-petrasvn-trunkunivie
Last change on this file since 02095e3 was e200b7a, checked in by Matthias Langer <matthias.langer@…>, 7 years ago

git-svn-id: http://flexpart.flexpart.eu:8088/svn/FlexPart90/trunk@4 ef8cc7e1-21b7-489e-abab-c1baa636049d

  • Property mode set to 100644
File size: 9.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  !                           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
996   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
234end subroutine interpol_wind
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG