Context Navigation

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
22	subroutine 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
145	end subroutine interpol_rain

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format