source: branches/jerome/src_flexwrf_v3.1/interpol_wind_short.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: 6.9 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_short(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)
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         *
40  !                                                                            *
41  !*****************************************************************************
42  !                                                                            *
43  ! Variables:                                                                 *
44  ! u,v,w              wind components                                         *
45  ! itime [s]          current temporal position                               *
46  ! memtime(3) [s]     times of the wind fields in memory                      *
47  ! xt,yt,zt           coordinates position for which wind data shall be       *
48  !                    calculated                                              *
49  !                                                                            *
50  ! Constants:                                                                 *
51  !                                                                            *
52  !*****************************************************************************
53
54  use par_mod
55  use com_mod
56!  use interpol_mod
57
58  implicit none
59
60  integer :: itime
61  real :: xt,yt,zt
62
63  ! Auxiliary variables needed for interpolation
64  real :: dz1,dz2,dz
65  real :: u1(2),v1(2),w1(2),uh(2),vh(2),wh(2)
66  integer :: i,m,n,indexh,indzh
67
68  real :: uprof(nzmax),vprof(nzmax),wprof(nzmax)
69  real :: usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax)
70  real :: rhoprof(nzmax),rhogradprof(nzmax)
71  real :: tkeprof(nzmax),pttprof(nzmax)
72  real :: u,v,w,usig,vsig,wsig,pvi
73
74  real :: p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2
75  integer :: ix,jy,ixp,jyp,ngrid,indz,indzp
76  logical :: depoindicator(maxspec)
77  logical :: indzindicator(nzmax)
78
79
80
81  !********************************************
82  ! Multilinear interpolation in time and space
83  !********************************************
84
85  ddx=xt-real(ix)
86  ddy=yt-real(jy)
87  rddx=1.-ddx
88  rddy=1.-ddy
89  p1=rddx*rddy
90  p2=ddx*rddy
91  p3=rddx*ddy
92  p4=ddx*ddy
93
94  ! Calculate variables for time interpolation
95  !*******************************************
96
97  dt1=real(itime-memtime(1))
98  dt2=real(memtime(2)-itime)
99  dtt=1./(dt1+dt2)
100
101  ! Determine the level below the current position for u,v
102  !*******************************************************
103
104  do i=2,nz
105    if (height(i).gt.zt) then
106      indz=i-1
107      goto 6
108    endif
109  end do
1106   continue
111
112
113  ! Vertical distance to the level below and above current position
114  !****************************************************************
115
116  dz=1./(height(indz+1)-height(indz))
117  dz1=(zt-height(indz))*dz
118  dz2=(height(indz+1)-zt)*dz
119
120
121  !**********************************************************************
122  ! 1.) Bilinear horizontal interpolation
123  ! This has to be done separately for 6 fields (Temporal(2)*Vertical(3))
124  !**********************************************************************
125
126  ! Loop over 2 time steps and 2 levels
127  !************************************
128
129  do m=1,2
130    indexh=memind(m)
131    do n=1,2
132      indzh=indz+n-1
133
134      if (ngrid.lt.0) then
135!!!$OMP CRITICAL
136        u1(n)=p1*uupol(ix ,jy ,indzh,indexh) &
137             +p2*uupol(ixp,jy ,indzh,indexh) &
138             +p3*uupol(ix ,jyp,indzh,indexh) &
139             +p4*uupol(ixp,jyp,indzh,indexh)
140        v1(n)=p1*vvpol(ix ,jy ,indzh,indexh) &
141             +p2*vvpol(ixp,jy ,indzh,indexh) &
142             +p3*vvpol(ix ,jyp,indzh,indexh) &
143             +p4*vvpol(ixp,jyp,indzh,indexh)
144!!!$OMP END CRITICAL
145      else
146!!!$OMP CRITICAL
147        u1(n)=p1*uu(ix ,jy ,indzh,indexh) &
148             +p2*uu(ixp,jy ,indzh,indexh) &
149             +p3*uu(ix ,jyp,indzh,indexh) &
150             +p4*uu(ixp,jyp,indzh,indexh)
151        v1(n)=p1*vv(ix ,jy ,indzh,indexh) &
152             +p2*vv(ixp,jy ,indzh,indexh) &
153             +p3*vv(ix ,jyp,indzh,indexh) &
154             +p4*vv(ixp,jyp,indzh,indexh)
155!!!$OMP END CRITICAL
156      endif
157!!!$OMP CRITICAL
158      w1(n)=p1*ww(ix ,jy ,indzh,indexh) &
159           +p2*ww(ixp,jy ,indzh,indexh) &
160           +p3*ww(ix ,jyp,indzh,indexh) &
161           +p4*ww(ixp,jyp,indzh,indexh)
162!!!$OMP END CRITICAL
163    end do
164
165
166  !**********************************
167  ! 2.) Linear vertical interpolation
168  !**********************************
169
170    uh(m)=dz2*u1(1)+dz1*u1(2)
171    vh(m)=dz2*v1(1)+dz1*v1(2)
172    wh(m)=dz2*w1(1)+dz1*w1(2)
173  end do
174
175
176
177  !************************************
178  ! 3.) Temporal interpolation (linear)
179  !************************************
180
181  u=(uh(1)*dt2+uh(2)*dt1)*dtt
182  v=(vh(1)*dt2+vh(2)*dt1)*dtt
183  w=(wh(1)*dt2+wh(2)*dt1)*dtt
184
185end subroutine interpol_wind_short
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG