source: trunk/src/interpol_wind_nests.f90

Last change on this file was 4, checked in by mlanger, 10 years ago
File size: 8.2 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_nests(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  !    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  ! u,v,w              wind components                                         *
42  ! itime [s]          current temporal position                               *
43  ! memtime(3) [s]     times of the wind fields in memory                      *
44  ! xt,yt,zt           coordinates position for which wind data shall be       *
45  !                    calculated                                              *
46  !                                                                            *
47  ! Constants:                                                                 *
48  !                                                                            *
49  !*****************************************************************************
50
51  use par_mod
52  use com_mod
53  use interpol_mod
54
55  implicit none
56
57  integer :: itime
58  real :: xt,yt,zt
59
60  ! Auxiliary variables needed for interpolation
61  real :: dz1,dz2,dz
62  real :: u1(2),v1(2),w1(2),uh(2),vh(2),wh(2)
63  real :: usl,vsl,wsl,usq,vsq,wsq,xaux
64  integer :: i,m,n,indexh,indzh
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  ! Determine the level below the current position for u,v
92  !*******************************************************
93
94  do i=2,nz
95    if (height(i).gt.zt) then
96      indz=i-1
97      goto 6
98    endif
99  end do
1006   continue
101
102
103  ! Vertical distance to the level below and above current position
104  !****************************************************************
105
106  dz=1./(height(indz+1)-height(indz))
107  dz1=(zt-height(indz))*dz
108  dz2=(height(indz+1)-zt)*dz
109
110
111  !**********************************************************************
112  ! 1.) Bilinear horizontal interpolation
113  ! This has to be done separately for 6 fields (Temporal(2)*Vertical(3))
114  !**********************************************************************
115
116  ! Loop over 2 time steps and 2 levels
117  !************************************
118
119  usl=0.
120  vsl=0.
121  wsl=0.
122  usq=0.
123  vsq=0.
124  wsq=0.
125  do m=1,2
126    indexh=memind(m)
127    do n=1,2
128      indzh=indz+n-1
129
130      u1(n)=p1*uun(ix ,jy ,indzh,indexh,ngrid) &
131           +p2*uun(ixp,jy ,indzh,indexh,ngrid) &
132           +p3*uun(ix ,jyp,indzh,indexh,ngrid) &
133           +p4*uun(ixp,jyp,indzh,indexh,ngrid)
134      v1(n)=p1*vvn(ix ,jy ,indzh,indexh,ngrid) &
135           +p2*vvn(ixp,jy ,indzh,indexh,ngrid) &
136           +p3*vvn(ix ,jyp,indzh,indexh,ngrid) &
137           +p4*vvn(ixp,jyp,indzh,indexh,ngrid)
138      w1(n)=p1*wwn(ix ,jy ,indzh,indexh,ngrid) &
139           +p2*wwn(ixp,jy ,indzh,indexh,ngrid) &
140           +p3*wwn(ix ,jyp,indzh,indexh,ngrid) &
141           +p4*wwn(ixp,jyp,indzh,indexh,ngrid)
142
143      usl=usl+uun(ix ,jy ,indzh,indexh,ngrid)+ &
144           uun(ixp,jy ,indzh,indexh,ngrid) &
145           +uun(ix ,jyp,indzh,indexh,ngrid)+ &
146           uun(ixp,jyp,indzh,indexh,ngrid)
147      vsl=vsl+vvn(ix ,jy ,indzh,indexh,ngrid)+ &
148           vvn(ixp,jy ,indzh,indexh,ngrid) &
149           +vvn(ix ,jyp,indzh,indexh,ngrid)+ &
150           vvn(ixp,jyp,indzh,indexh,ngrid)
151      wsl=wsl+wwn(ix ,jy ,indzh,indexh,ngrid)+ &
152           wwn(ixp,jy ,indzh,indexh,ngrid) &
153           +wwn(ix ,jyp,indzh,indexh,ngrid)+ &
154           wwn(ixp,jyp,indzh,indexh,ngrid)
155
156      usq=usq+uun(ix ,jy ,indzh,indexh,ngrid)* &
157           uun(ix ,jy ,indzh,indexh,ngrid)+ &
158           uun(ixp,jy ,indzh,indexh,ngrid)*uun(ixp,jy ,indzh,indexh,ngrid)+ &
159           uun(ix ,jyp,indzh,indexh,ngrid)*uun(ix ,jyp,indzh,indexh,ngrid)+ &
160           uun(ixp,jyp,indzh,indexh,ngrid)*uun(ixp,jyp,indzh,indexh,ngrid)
161      vsq=vsq+vvn(ix ,jy ,indzh,indexh,ngrid)* &
162           vvn(ix ,jy ,indzh,indexh,ngrid)+ &
163           vvn(ixp,jy ,indzh,indexh,ngrid)*vvn(ixp,jy ,indzh,indexh,ngrid)+ &
164           vvn(ix ,jyp,indzh,indexh,ngrid)*vvn(ix ,jyp,indzh,indexh,ngrid)+ &
165           vvn(ixp,jyp,indzh,indexh,ngrid)*vvn(ixp,jyp,indzh,indexh,ngrid)
166      wsq=wsq+wwn(ix ,jy ,indzh,indexh,ngrid)* &
167           wwn(ix ,jy ,indzh,indexh,ngrid)+ &
168           wwn(ixp,jy ,indzh,indexh,ngrid)*wwn(ixp,jy ,indzh,indexh,ngrid)+ &
169           wwn(ix ,jyp,indzh,indexh,ngrid)*wwn(ix ,jyp,indzh,indexh,ngrid)+ &
170           wwn(ixp,jyp,indzh,indexh,ngrid)*wwn(ixp,jyp,indzh,indexh,ngrid)
171    end do
172
173
174  !**********************************
175  ! 2.) Linear vertical interpolation
176  !**********************************
177
178    uh(m)=dz2*u1(1)+dz1*u1(2)
179    vh(m)=dz2*v1(1)+dz1*v1(2)
180    wh(m)=dz2*w1(1)+dz1*w1(2)
181  end do
182
183
184  !************************************
185  ! 3.) Temporal interpolation (linear)
186  !************************************
187
188  u=(uh(1)*dt2+uh(2)*dt1)*dtt
189  v=(vh(1)*dt2+vh(2)*dt1)*dtt
190  w=(wh(1)*dt2+wh(2)*dt1)*dtt
191
192
193  ! Compute standard deviations
194  !****************************
195
196  xaux=usq-usl*usl/16.
197  if (xaux.lt.eps) then
198    usig=0.
199  else
200    usig=sqrt(xaux/15.)
201  endif
202
203  xaux=vsq-vsl*vsl/16.
204  if (xaux.lt.eps) then
205    vsig=0.
206  else
207    vsig=sqrt(xaux/15.)
208  endif
209
210
211  xaux=wsq-wsl*wsl/16.
212  if (xaux.lt.eps) then
213    wsig=0.
214  else
215    wsig=sqrt(xaux/15.)
216  endif
217
218end subroutine interpol_wind_nests
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG