CTBTO FLEXPART WO4 (2015-10-15)
 All Classes Files Functions Variables
interpol_rain_nests.f90
Go to the documentation of this file.
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 
22 subroutine interpol_rain_nests(yy1,yy2,yy3,nxmaxn,nymaxn,nzmax, &
23  maxnests,ngrid,nxn,nyn,memind,xt,yt,level,itime1,itime2,itime, &
24  yint1,yint2,yint3)
25  ! i i i i i i
26  ! i i i i i i i i i i i
27  ! o o o
28  !****************************************************************************
29  ! *
30  ! Interpolation of meteorological fields on 2-d model layers for nested *
31  ! grids. This routine is related to levlin3interpol.f for the mother domain*
32  ! *
33  ! In horizontal direction bilinear interpolation interpolation is used. *
34  ! Temporally a linear interpolation is used. *
35  ! Three fields are interpolated at the same time. *
36  ! *
37  ! This is a special version of levlininterpol to save CPU time. *
38  ! *
39  ! 1 first time *
40  ! 2 second time *
41  ! *
42  ! *
43  ! Author: A. Stohl *
44  ! *
45  ! 15 March 2000 *
46  ! *
47  !****************************************************************************
48  ! *
49  ! Variables: *
50  ! *
51  ! dt1,dt2 time differences between fields and current position *
52  ! dz1,dz2 z distance between levels and current position *
53  ! height(nzmax) heights of the model levels *
54  ! indexh help variable *
55  ! indz the level closest to the current trajectory position *
56  ! indzh help variable *
57  ! itime current time *
58  ! itime1 time of the first wind field *
59  ! itime2 time of the second wind field *
60  ! ix,jy x,y coordinates of lower left subgrid point *
61  ! level level at which interpolation shall be done *
62  ! memind(3) points to the places of the wind fields *
63  ! nx,ny actual field dimensions in x,y and z direction *
64  ! nxmax,nymax,nzmax maximum field dimensions in x,y and z direction *
65  ! xt current x coordinate *
66  ! yint the final interpolated value *
67  ! yt current y coordinate *
68  ! yy(0:nxmax,0:nymax,nzmax,3) meteorological field used for interpolation *
69  ! zt current z coordinate *
70  ! *
71  !****************************************************************************
72 
73  implicit none
74 
75  integer :: maxnests,ngrid
76  integer :: nxn(maxnests),nyn(maxnests),nxmaxn,nymaxn,nzmax,memind(2)
77  integer :: m,ix,jy,ixp,jyp,itime,itime1,itime2,level,indexh
78  real :: yy1(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
79  real :: yy2(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
80  real :: yy3(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
81  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
82  real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
83 
84 
85 
86  ! If point at border of grid -> small displacement into grid
87  !***********************************************************
88 
89  if (xt.ge.(real(nxn(ngrid)-1)-0.0001)) &
90  xt=real(nxn(ngrid)-1)-0.0001
91  if (yt.ge.(real(nyn(ngrid)-1)-0.0001)) &
92  yt=real(nyn(ngrid)-1)-0.0001
93 
94 
95 
96  !**********************************************************************
97  ! 1.) Bilinear horizontal interpolation
98  ! This has to be done separately for 2 fields (Temporal)
99  !*******************************************************
100 
101  ! Determine the lower left corner and its distance to the current position
102  !*************************************************************************
103 
104  ix=int(xt)
105  jy=int(yt)
106  ixp=ix+1
107  jyp=jy+1
108  ddx=xt-real(ix)
109  ddy=yt-real(jy)
110  rddx=1.-ddx
111  rddy=1.-ddy
112  p1=rddx*rddy
113  p2=ddx*rddy
114  p3=rddx*ddy
115  p4=ddx*ddy
116 
117 
118  ! Loop over 2 time steps
119  !***********************
120 
121  do m=1,2
122  indexh=memind(m)
123 
124  y1(m)=p1*yy1(ix ,jy ,level,indexh,ngrid) &
125  + p2*yy1(ixp,jy ,level,indexh,ngrid) &
126  + p3*yy1(ix ,jyp,level,indexh,ngrid) &
127  + p4*yy1(ixp,jyp,level,indexh,ngrid)
128  y2(m)=p1*yy2(ix ,jy ,level,indexh,ngrid) &
129  + p2*yy2(ixp,jy ,level,indexh,ngrid) &
130  + p3*yy2(ix ,jyp,level,indexh,ngrid) &
131  + p4*yy2(ixp,jyp,level,indexh,ngrid)
132  y3(m)=p1*yy3(ix ,jy ,level,indexh,ngrid) &
133  + p2*yy3(ixp,jy ,level,indexh,ngrid) &
134  + p3*yy3(ix ,jyp,level,indexh,ngrid) &
135  + p4*yy3(ixp,jyp,level,indexh,ngrid)
136  end do
137 
138 
139  !************************************
140  ! 2.) Temporal interpolation (linear)
141  !************************************
142 
143  dt1=real(itime-itime1)
144  dt2=real(itime2-itime)
145  dt=dt1+dt2
146 
147  yint1=(y1(1)*dt2+y1(2)*dt1)/dt
148  yint2=(y2(1)*dt2+y2(2)*dt1)/dt
149  yint3=(y3(1)*dt2+y3(2)*dt1)/dt
150 
151 
152 end subroutine interpol_rain_nests
subroutine interpol_rain_nests(yy1, yy2, yy3, nxmaxn, nymaxn, nzmax, maxnests, ngrid, nxn, nyn, memind, xt, yt, level, itime1, itime2, itime, yint1, yint2, yint3)