source: branches/jerome/src_flexwrf_v3.1/interpol_all.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 10 years ago

sources for flexwrf v3.1

File size: 12.2 KB
Line 
1!***********************************************************************
2!* Copyright 2012,2013                                                *
3!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine,       *
4!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso,          *
6!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann,   *
7!*                                                                     *
8!* This file is part of FLEXPART WRF                                   *
9!*                                                                     *
10!* FLEXPART is free software: you can redistribute it and/or modify    *
11!* it under the terms of the GNU General Public License as published by*
12!* the Free Software Foundation, either version 3 of the License, or   *
13!* (at your option) any later version.                                 *
14!*                                                                     *
15!* FLEXPART is distributed in the hope that it will be useful,         *
16!* but WITHOUT ANY WARRANTY; without even the implied warranty of      *
17!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
18!* GNU General Public License for more details.                        *
19!*                                                                     *
20!* You should have received a copy of the GNU General Public License   *
21!* along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
22!***********************************************************************
23      subroutine interpol_all(itime,xt,yt,zt, &
24    uprof,vprof,wprof, usigprof,vsigprof,wsigprof, &
25    rhoprof,rhogradprof, tkeprof,pttprof, &
26    u,v,w,usig,vsig,wsig,pvi, &
27    p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2, &
28    ix,jy,ixp,jyp,ngrid,indz,indzp, depoindicator, &
29    indzindicator, &
30   ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
31   sigw,dsigwdz,dsigw2dz,mu,mv)
32
33!                               i   i  i  i
34!*******************************************************************************
35!                                                                              *
36!  This subroutine interpolates everything that is needed for calculating the  *
37!  dispersion.                                                                 *
38!                                                                              *
39!    Author: A. Stohl                                                          *
40!                                                                              *
41!    16 December 1997                                                          *
42!                                                                              *
43!  Revision March 2005 by AST : all output variables in common block           *
44!                               calculation of standard deviation done in this *
45!                               routine rather than subroutine call in order   *
46!                               to save computation time                       *
47!                                                                              *
48!*******************************************************************************
49!                                                                              *
50! Variables:                                                                   *
51! itime [s]          current temporal position                                 *
52! memtime(3) [s]     times of the wind fields in memory                        *
53! xt,yt,zt           coordinates position for which wind data shall be calculat*
54!                                                                              *
55! Constants:                                                                   *
56!                                                                              *
57!*******************************************************************************
58!     12 JUNE 2007, compute tkeprof,   y4
59!     25 June 2007, compute pttprof,   y5
60!                   compute tkeprof for all levels
61!*******************************************************************************
62  use par_mod
63  use com_mod
64!  use interpol_mod
65!  use hanna_mod
66  implicit none
67
68      integer :: itime
69      real :: xt,yt,zt
70
71! Auxiliary variables needed for interpolation
72      real :: ust1(2),wst1(2),oli1(2),oliaux
73      real :: y1(2),y2(2),y3(2),rho1(2),rhograd1(2),y4(2),y5(2)
74      real :: usl,vsl,wsl,usq,vsq,wsq,xaux
75      integer :: i,m,n,indexh,n2
76      real :: tkeprof2
77
78      real,parameter ::eps=1.0e-30
79
80  real :: uprof(nzmax),vprof(nzmax),wprof(nzmax)
81  real :: usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax)
82  real :: rhoprof(nzmax),rhogradprof(nzmax)
83  real :: tkeprof(nzmax),pttprof(nzmax)
84  real :: u,v,w,usig,vsig,wsig,pvi
85
86  real :: p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2
87  integer :: ix,jy,ixp,jyp,ngrid,indz,indzp
88  logical :: depoindicator(maxspec)
89  logical :: indzindicator(nzmax)
90  real :: ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw
91  real :: sigw,dsigwdz,dsigw2dz,mu,mv
92
93
94!********************************************
95! Multilinear interpolation in time and space
96!********************************************
97
98! Determine the lower left corner and its distance to the current position
99!*************************************************************************
100
101      ddx=xt-real(ix)                     
102      ddy=yt-real(jy)
103      rddx=1.-ddx
104      rddy=1.-ddy
105      p1=rddx*rddy
106      p2=ddx*rddy
107      p3=rddx*ddy
108      p4=ddx*ddy
109
110! Calculate variables for time interpolation
111!*******************************************
112
113      dt1=real(itime-memtime(1))
114      dt2=real(memtime(2)-itime)
115      dtt=1./(dt1+dt2)
116
117
118!*****************************************
119! 1. Interpolate u*, w* and Obukhov length
120!*****************************************
121
122! a) Bilinear horizontal interpolation
123
124  do m=1,2
125    indexh=memind(m)
126
127    ust1(m)=p1*ustar(ix ,jy ,1,indexh) &
128         + p2*ustar(ixp,jy ,1,indexh) &
129         + p3*ustar(ix ,jyp,1,indexh) &
130         + p4*ustar(ixp,jyp,1,indexh)
131    wst1(m)=p1*wstar(ix ,jy ,1,indexh) &
132         + p2*wstar(ixp,jy ,1,indexh) &
133         + p3*wstar(ix ,jyp,1,indexh) &
134         + p4*wstar(ixp,jyp,1,indexh)
135    oli1(m)=p1*oli(ix ,jy ,1,indexh) &
136         + p2*oli(ixp,jy ,1,indexh) &
137         + p3*oli(ix ,jyp,1,indexh) &
138         + p4*oli(ixp,jyp,1,indexh)
139  end do
140       mu =p1*m_x(ix ,jy ,1) &
141         + p2*m_x(ixp,jy ,1) &
142         + p3*m_x(ix ,jyp,1) &
143         + p4*m_x(ixp,jyp,1)
144       mv =p1*m_y(ix ,jy ,1) &
145         + p2*m_y(ixp,jy ,1) &
146         + p3*m_y(ix ,jyp,1) &
147         + p4*m_y(ixp,jyp,1)
148
149! b) Temporal interpolation
150
151      ust=(ust1(1)*dt2+ust1(2)*dt1)*dtt
152      wst=(wst1(1)*dt2+wst1(2)*dt1)*dtt
153      oliaux=(oli1(1)*dt2+oli1(2)*dt1)*dtt
154
155      if (oliaux.ne.0.) then
156        ol=1./oliaux
157      else
158        ol=99999.
159      endif
160
161
162!*****************************************************
163! 2. Interpolate vertical profiles of u,v,w,rho,drhodz
164!*****************************************************
165
166
167! Determine the level below the current position
168!***********************************************
169
170  do i=2,nz
171    if (height(i).gt.zt) then
172      indz=i-1
173      indzp=i
174      goto 6
175    endif
176  end do
1776   continue
178
179!**************************************
180! 1.) Bilinear horizontal interpolation
181! 2.) Temporal interpolation (linear)
182!**************************************
183
184! Loop over 2 time steps and indz levels
185!***************************************
186
187  do n=indz,indzp
188    usl=0.
189    vsl=0.
190    wsl=0.
191    usq=0.
192    vsq=0.
193    wsq=0.
194    do m=1,2
195      indexh=memind(m)
196      if (ngrid.lt.0) then
197        y1(m)=p1*uupol(ix ,jy ,n,indexh) &
198             +p2*uupol(ixp,jy ,n,indexh) &
199             +p3*uupol(ix ,jyp,n,indexh) &
200             +p4*uupol(ixp,jyp,n,indexh)
201        y2(m)=p1*vvpol(ix ,jy ,n,indexh) &
202             +p2*vvpol(ixp,jy ,n,indexh) &
203             +p3*vvpol(ix ,jyp,n,indexh) &
204             +p4*vvpol(ixp,jyp,n,indexh)
205        usl=usl+uupol(ix ,jy ,n,indexh)+uupol(ixp,jy ,n,indexh) &
206             +uupol(ix ,jyp,n,indexh)+uupol(ixp,jyp,n,indexh)
207        vsl=vsl+vvpol(ix ,jy ,n,indexh)+vvpol(ixp,jy ,n,indexh) &
208             +vvpol(ix ,jyp,n,indexh)+vvpol(ixp,jyp,n,indexh)
209
210        usq=usq+uupol(ix ,jy ,n,indexh)*uupol(ix ,jy ,n,indexh)+ &
211             uupol(ixp,jy ,n,indexh)*uupol(ixp,jy ,n,indexh)+ &
212             uupol(ix ,jyp,n,indexh)*uupol(ix ,jyp,n,indexh)+ &
213             uupol(ixp,jyp,n,indexh)*uupol(ixp,jyp,n,indexh)
214        vsq=vsq+vvpol(ix ,jy ,n,indexh)*vvpol(ix ,jy ,n,indexh)+ &
215             vvpol(ixp,jy ,n,indexh)*vvpol(ixp,jy ,n,indexh)+ &
216             vvpol(ix ,jyp,n,indexh)*vvpol(ix ,jyp,n,indexh)+ &
217             vvpol(ixp,jyp,n,indexh)*vvpol(ixp,jyp,n,indexh)
218          else
219        y1(m)=p1*uu(ix ,jy ,n,indexh) &
220             +p2*uu(ixp,jy ,n,indexh) &
221             +p3*uu(ix ,jyp,n,indexh) &
222             +p4*uu(ixp,jyp,n,indexh)
223        y2(m)=p1*vv(ix ,jy ,n,indexh) &
224             +p2*vv(ixp,jy ,n,indexh) &
225             +p3*vv(ix ,jyp,n,indexh) &
226             +p4*vv(ixp,jyp,n,indexh)
227        usl=usl+uu(ix ,jy ,n,indexh)+uu(ixp,jy ,n,indexh) &
228             +uu(ix ,jyp,n,indexh)+uu(ixp,jyp,n,indexh)
229        vsl=vsl+vv(ix ,jy ,n,indexh)+vv(ixp,jy ,n,indexh) &
230             +vv(ix ,jyp,n,indexh)+vv(ixp,jyp,n,indexh)
231
232        usq=usq+uu(ix ,jy ,n,indexh)*uu(ix ,jy ,n,indexh)+ &
233             uu(ixp,jy ,n,indexh)*uu(ixp,jy ,n,indexh)+ &
234             uu(ix ,jyp,n,indexh)*uu(ix ,jyp,n,indexh)+ &
235             uu(ixp,jyp,n,indexh)*uu(ixp,jyp,n,indexh)
236        vsq=vsq+vv(ix ,jy ,n,indexh)*vv(ix ,jy ,n,indexh)+ &
237             vv(ixp,jy ,n,indexh)*vv(ixp,jy ,n,indexh)+ &
238             vv(ix ,jyp,n,indexh)*vv(ix ,jyp,n,indexh)+ &
239             vv(ixp,jyp,n,indexh)*vv(ixp,jyp,n,indexh)
240          endif
241      y3(m)=p1*ww(ix ,jy ,n,indexh) &
242           +p2*ww(ixp,jy ,n,indexh) &
243           +p3*ww(ix ,jyp,n,indexh) &
244           +p4*ww(ixp,jyp,n,indexh)
245      rhograd1(m)=p1*drhodz(ix ,jy ,n,indexh) &
246           +p2*drhodz(ixp,jy ,n,indexh) &
247           +p3*drhodz(ix ,jyp,n,indexh) &
248           +p4*drhodz(ixp,jyp,n,indexh)
249      rho1(m)=p1*rho(ix ,jy ,n,indexh) &
250           +p2*rho(ixp,jy ,n,indexh) &
251           +p3*rho(ix ,jyp,n,indexh) &
252           +p4*rho(ixp,jyp,n,indexh)
253      wsl=wsl+ww(ix ,jy ,n,indexh)+ww(ixp,jy ,n,indexh) &
254           +ww(ix ,jyp,n,indexh)+ww(ixp,jyp,n,indexh)
255      wsq=wsq+ww(ix ,jy ,n,indexh)*ww(ix ,jy ,n,indexh)+ &
256           ww(ixp,jy ,n,indexh)*ww(ixp,jy ,n,indexh)+ &
257           ww(ix ,jyp,n,indexh)*ww(ix ,jyp,n,indexh)+ &
258           ww(ixp,jyp,n,indexh)*ww(ixp,jyp,n,indexh)
259
260
261         enddo
262        uprof(n)=(y1(1)*dt2+y1(2)*dt1)*dtt
263        vprof(n)=(y2(1)*dt2+y2(2)*dt1)*dtt
264        wprof(n)=(y3(1)*dt2+y3(2)*dt1)*dtt
265        rhoprof(n)=(rho1(1)*dt2+rho1(2)*dt1)*dtt
266        rhogradprof(n)=(rhograd1(1)*dt2+rhograd1(2)*dt1)*dtt
267        indzindicator(n)=.false.
268
269! Compute standard deviations
270!****************************
271
272        xaux=usq-usl*usl/8.
273        if (xaux.lt.eps) then
274          usigprof(n)=0.
275        else
276          usigprof(n)=sqrt(xaux/7.)
277        endif
278
279        xaux=vsq-vsl*vsl/8.
280        if (xaux.lt.eps) then
281          vsigprof(n)=0.
282        else
283          vsigprof(n)=sqrt(xaux/7.)
284        endif
285
286
287        xaux=wsq-wsl*wsl/8.
288        if (xaux.lt.eps) then
289          wsigprof(n)=0.
290        else
291          wsigprof(n)=sqrt(xaux/7.)
292        endif
293
294        enddo
295
296! compute TKE for all levels, used for estimating turb length scale
297         tkeprof2=-999.
298!        tkeprof2=0.
299        do n=1,nz
300         do m=1,2
301          indexh=memind(m)
302         y4(m) =p1*tke(ix ,jy ,n,indexh) &
303                 +p2*tke(ixp,jy ,n,indexh) &
304                 +p3*tke(ix ,jyp,n,indexh) &
305                 +p4*tke(ixp,jyp,n,indexh)
306          y5(m) =p1*ptt(ix ,jy ,n,indexh) &
307                 +p2*ptt(ixp,jy ,n,indexh) &
308                 +p3*ptt(ix ,jyp,n,indexh) &
309                 +p4*ptt(ixp,jyp,n,indexh)
310          enddo
311       tkeprof(n)=(y4(1)*dt2+y4(2)*dt1)*dtt
312       pttprof(n)=(y5(1)*dt2+y5(2)*dt1)*dtt
313!        if (tkeprof(n).gt.tkeprof2) then
314!         tkeprof2=tkeprof(n)
315!         n2=n
316!         endif
317!         if (n.lt.20) then
318!         tkeprof2=tkeprof2+tkeprof(n)
319!         n2=20
320!         endif       
321        enddo
322!           tkeprof(1)=0.33*tkeprof(1)+0.33*tkeprof(2)+0.33*tkeprof(3)
323!         write(*,*)'interpol_all,itime,xt,yt,zt',itime,xt,yt,zt
324!         write(*,*)(tkeprof(n),n=1,nz)
325
326end subroutine interpol_all
327
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG