source: trunk/src/interpol_rain.f90 @ 24

Last change on this file since 24 was 24, checked in by igpis, 10 years ago

version 9.2 beta. Changes from HH, AST, MC, NIK, IP. Changes in vert transform. New SPECIES input includes scavenging coefficients

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