source: flexpart.git/src/interpol_wind_nests.f90

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file was 92fab65, checked in by Ignacio Pisso <ip@…>, 4 years ago

add SPDX-License-Identifier to all .f90 files

  • Property mode set to 100644
File size: 6.9 KB
Line 
1! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
2! SPDX-License-Identifier: GPL-3.0-or-later
3
4subroutine interpol_wind_nests(itime,xt,yt,zt)
5  !                                 i   i  i  i
6  !*****************************************************************************
7  !                                                                            *
8  !  This subroutine interpolates the wind data to current trajectory position.*
9  !                                                                            *
10  !    Author: A. Stohl                                                        *
11  !                                                                            *
12  !    16 December 1997                                                        *
13  !    16 December 1997                                                        *
14  !                                                                            *
15  !  Revision March 2005 by AST : all output variables in common block cal-    *
16  !                               culation of standard deviation done in this  *
17  !                               routine rather than subroutine call in order *
18  !                               to save computation time                     *
19  !                                                                            *
20  !*****************************************************************************
21  !                                                                            *
22  ! Variables:                                                                 *
23  ! u,v,w              wind components                                         *
24  ! itime [s]          current temporal position                               *
25  ! memtime(3) [s]     times of the wind fields in memory                      *
26  ! xt,yt,zt           coordinates position for which wind data shall be       *
27  !                    calculated                                              *
28  !                                                                            *
29  ! Constants:                                                                 *
30  !                                                                            *
31  !*****************************************************************************
32
33  use par_mod
34  use com_mod
35  use interpol_mod
36
37  implicit none
38
39  integer :: itime
40  real :: xt,yt,zt
41
42  ! Auxiliary variables needed for interpolation
43  real :: dz1,dz2,dz
44  real :: u1(2),v1(2),w1(2),uh(2),vh(2),wh(2)
45  real :: usl,vsl,wsl,usq,vsq,wsq,xaux
46  integer :: i,m,n,indexh,indzh
47  real,parameter :: eps=1.0e-30
48
49
50  !********************************************
51  ! Multilinear interpolation in time and space
52  !********************************************
53
54  ! Determine the lower left corner and its distance to the current position
55  !*************************************************************************
56
57  ddx=xt-real(ix)
58  ddy=yt-real(jy)
59  rddx=1.-ddx
60  rddy=1.-ddy
61  p1=rddx*rddy
62  p2=ddx*rddy
63  p3=rddx*ddy
64  p4=ddx*ddy
65
66  ! Calculate variables for time interpolation
67  !*******************************************
68
69  dt1=real(itime-memtime(1))
70  dt2=real(memtime(2)-itime)
71  dtt=1./(dt1+dt2)
72
73  ! Determine the level below the current position for u,v
74  !*******************************************************
75
76  do i=2,nz
77    if (height(i).gt.zt) then
78      indz=i-1
79      goto 6
80    endif
81  end do
826   continue
83
84
85  ! Vertical distance to the level below and above current position
86  !****************************************************************
87
88  dz=1./(height(indz+1)-height(indz))
89  dz1=(zt-height(indz))*dz
90  dz2=(height(indz+1)-zt)*dz
91
92
93  !**********************************************************************
94  ! 1.) Bilinear horizontal interpolation
95  ! This has to be done separately for 6 fields (Temporal(2)*Vertical(3))
96  !**********************************************************************
97
98  ! Loop over 2 time steps and 2 levels
99  !************************************
100
101  usl=0.
102  vsl=0.
103  wsl=0.
104  usq=0.
105  vsq=0.
106  wsq=0.
107  do m=1,2
108    indexh=memind(m)
109    do n=1,2
110      indzh=indz+n-1
111
112      u1(n)=p1*uun(ix ,jy ,indzh,indexh,ngrid) &
113           +p2*uun(ixp,jy ,indzh,indexh,ngrid) &
114           +p3*uun(ix ,jyp,indzh,indexh,ngrid) &
115           +p4*uun(ixp,jyp,indzh,indexh,ngrid)
116      v1(n)=p1*vvn(ix ,jy ,indzh,indexh,ngrid) &
117           +p2*vvn(ixp,jy ,indzh,indexh,ngrid) &
118           +p3*vvn(ix ,jyp,indzh,indexh,ngrid) &
119           +p4*vvn(ixp,jyp,indzh,indexh,ngrid)
120      w1(n)=p1*wwn(ix ,jy ,indzh,indexh,ngrid) &
121           +p2*wwn(ixp,jy ,indzh,indexh,ngrid) &
122           +p3*wwn(ix ,jyp,indzh,indexh,ngrid) &
123           +p4*wwn(ixp,jyp,indzh,indexh,ngrid)
124
125      usl=usl+uun(ix ,jy ,indzh,indexh,ngrid)+ &
126           uun(ixp,jy ,indzh,indexh,ngrid) &
127           +uun(ix ,jyp,indzh,indexh,ngrid)+ &
128           uun(ixp,jyp,indzh,indexh,ngrid)
129      vsl=vsl+vvn(ix ,jy ,indzh,indexh,ngrid)+ &
130           vvn(ixp,jy ,indzh,indexh,ngrid) &
131           +vvn(ix ,jyp,indzh,indexh,ngrid)+ &
132           vvn(ixp,jyp,indzh,indexh,ngrid)
133      wsl=wsl+wwn(ix ,jy ,indzh,indexh,ngrid)+ &
134           wwn(ixp,jy ,indzh,indexh,ngrid) &
135           +wwn(ix ,jyp,indzh,indexh,ngrid)+ &
136           wwn(ixp,jyp,indzh,indexh,ngrid)
137
138      usq=usq+uun(ix ,jy ,indzh,indexh,ngrid)* &
139           uun(ix ,jy ,indzh,indexh,ngrid)+ &
140           uun(ixp,jy ,indzh,indexh,ngrid)*uun(ixp,jy ,indzh,indexh,ngrid)+ &
141           uun(ix ,jyp,indzh,indexh,ngrid)*uun(ix ,jyp,indzh,indexh,ngrid)+ &
142           uun(ixp,jyp,indzh,indexh,ngrid)*uun(ixp,jyp,indzh,indexh,ngrid)
143      vsq=vsq+vvn(ix ,jy ,indzh,indexh,ngrid)* &
144           vvn(ix ,jy ,indzh,indexh,ngrid)+ &
145           vvn(ixp,jy ,indzh,indexh,ngrid)*vvn(ixp,jy ,indzh,indexh,ngrid)+ &
146           vvn(ix ,jyp,indzh,indexh,ngrid)*vvn(ix ,jyp,indzh,indexh,ngrid)+ &
147           vvn(ixp,jyp,indzh,indexh,ngrid)*vvn(ixp,jyp,indzh,indexh,ngrid)
148      wsq=wsq+wwn(ix ,jy ,indzh,indexh,ngrid)* &
149           wwn(ix ,jy ,indzh,indexh,ngrid)+ &
150           wwn(ixp,jy ,indzh,indexh,ngrid)*wwn(ixp,jy ,indzh,indexh,ngrid)+ &
151           wwn(ix ,jyp,indzh,indexh,ngrid)*wwn(ix ,jyp,indzh,indexh,ngrid)+ &
152           wwn(ixp,jyp,indzh,indexh,ngrid)*wwn(ixp,jyp,indzh,indexh,ngrid)
153    end do
154
155
156  !**********************************
157  ! 2.) Linear vertical interpolation
158  !**********************************
159
160    uh(m)=dz2*u1(1)+dz1*u1(2)
161    vh(m)=dz2*v1(1)+dz1*v1(2)
162    wh(m)=dz2*w1(1)+dz1*w1(2)
163  end do
164
165
166  !************************************
167  ! 3.) Temporal interpolation (linear)
168  !************************************
169
170  u=(uh(1)*dt2+uh(2)*dt1)*dtt
171  v=(vh(1)*dt2+vh(2)*dt1)*dtt
172  w=(wh(1)*dt2+wh(2)*dt1)*dtt
173
174
175  ! Compute standard deviations
176  !****************************
177
178  xaux=usq-usl*usl/16.
179  if (xaux.lt.eps) then
180    usig=0.
181  else
182    usig=sqrt(xaux/15.)
183  endif
184
185  xaux=vsq-vsl*vsl/16.
186  if (xaux.lt.eps) then
187    vsig=0.
188  else
189    vsig=sqrt(xaux/15.)
190  endif
191
192
193  xaux=wsq-wsl*wsl/16.
194  if (xaux.lt.eps) then
195    wsig=0.
196  else
197    wsig=sqrt(xaux/15.)
198  endif
199
200end subroutine interpol_wind_nests
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG