Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

interpol_rain.f90 @ 14

Last change on this file since 14 was 14, checked in by igpis, 11 years ago
based on 9.1from hasod; 9.1.1 add sack (trunk); 9.1.2 add NIK scavenging; 9.1.3 add new pesei depo scheme; 9.1.4 warning at readpositions; 9.1.5 xuekun FNL; 9.1.6 rlt FLEXINVERT; 9.1.7 update dates, check slash in COMMAND
File size: 9.4 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	subroutine interpol_rain(yy1,yy2,yy3,iy1,iy2,nxmax,nymax,nzmax,nx, &
28	ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, &
29	intiy1,intiy2,icmv)
30	! i i i i i i i
31	! i i i i i i i i o o o
32
33	!****************************************************************************
34	! *
35	! Interpolation of meteorological fields on 2-d model layers. *
36	! In horizontal direction bilinear interpolation interpolation is used. *
37	! Temporally a linear interpolation is used. *
38	! Three fields are interpolated at the same time. *
39	! *
40	! This is a special version of levlininterpol to save CPU time. *
41	! *
42	! 1 first time *
43	! 2 second time *
44	! *
45	! *
46	! Author: A. Stohl *
47	! *
48	! 30 August 1996 *
49	!
50	!* Petra Seibert, 2011/2012:
51	!* Add interpolation of cloud bottom and cloud thickness
52	!* for fix to SE's new wet scavenging scheme *
53	!****************************************************************************
54	! *
55	! Variables: *
56	! *
57	! dt1,dt2 time differences between fields and current position *
58	! dz1,dz2 z distance between levels and current position *
59	! height(nzmax) heights of the model levels *
60	! indexh help variable *
61	! indz the level closest to the current trajectory position *
62	! indzh help variable *
63	! itime current time *
64	! itime1 time of the first wind field *
65	! itime2 time of the second wind field *
66	! ix,jy x,y coordinates of lower left subgrid point *
67	! level level at which interpolation shall be done *
68	! memind(3) points to the places of the wind fields *
69	! nx,ny actual field dimensions in x,y and z direction *
70	! nxmax,nymax,nzmax maximum field dimensions in x,y and z direction *
71	! xt current x coordinate *
72	! yint the final interpolated value *
73	! yt current y coordinate *
74	! yy(0:nxmax,0:nymax,nzmax,3) meteorological field used for interpolation *
75	! zt current z coordinate *
76	! *
77	!****************************************************************************
78
79	implicit none
80
81	integer :: nx,ny,nxmax,nymax,nzmax,memind(2),m,ix,jy,ixp,jyp
82	!integer :: itime,itime1,itime2,level,indexh
83	integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4
84	integer :: intiy1,intiy2,ipsum,icmv
85	real :: yy1(0:nxmax-1,0:nymax-1,nzmax,2)
86	real :: yy2(0:nxmax-1,0:nymax-1,nzmax,2)
87	real :: yy3(0:nxmax-1,0:nymax-1,nzmax,2)
88	integer iy1(0:nxmax-1,0:nymax-1,2),iy2(0:nxmax-1,0:nymax-1,2)
89	real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),yi1(2),yi2(2)
90	!real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
91	real :: xt,yt,yint1,yint2,yint3,yint4,p1,p2,p3,p4
92	!real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
93
94
95
96	! If point at border of grid -> small displacement into grid
97	!***********************************************************
98
99	if (xt.ge.real(nx-1)) xt=real(nx-1)-0.00001
100	if (yt.ge.real(ny-1)) yt=real(ny-1)-0.00001
101
102
103
104	!**********************************************************************
105	! 1.) Bilinear horizontal interpolation
106	! This has to be done separately for 2 fields (Temporal)
107	!*******************************************************
108
109	! Determine the lower left corner and its distance to the current position
110	!*************************************************************************
111
112	ix=int(xt)
113	jy=int(yt)
114	ixp=ix+1
115	jyp=jy+1
116	ddx=xt-real(ix)
117	ddy=yt-real(jy)
118	rddx=1.-ddx
119	rddy=1.-ddy
120	p1=rddx*rddy
121	p2=ddx*rddy
122	p3=rddx*ddy
123	p4=ddx*ddy
124
125
126	! Loop over 2 time steps
127	!***********************
128
129	do m=1,2
130	indexh=memind(m)
131
132	y1(m)=p1*yy1(ix ,jy ,level,indexh) &
133	+ p2*yy1(ixp,jy ,level,indexh) &
134	+ p3*yy1(ix ,jyp,level,indexh) &
135	+ p4*yy1(ixp,jyp,level,indexh)
136	y2(m)=p1*yy2(ix ,jy ,level,indexh) &
137	+ p2*yy2(ixp,jy ,level,indexh) &
138	+ p3*yy2(ix ,jyp,level,indexh) &
139	+ p4*yy2(ixp,jyp,level,indexh)
140	y3(m)=p1*yy3(ix ,jy ,level,indexh) &
141	+ p2*yy3(ixp,jy ,level,indexh) &
142	+ p3*yy3(ix ,jyp,level,indexh) &
143	+ p4*yy3(ixp,jyp,level,indexh)
144
145	!CPS clouds:
146	ip1=1
147	ip2=1
148	ip3=1
149	ip4=1
150	if (iy1(ix ,jy ,indexh) .eq. icmv) ip1=0
151	if (iy1(ixp,jy ,indexh) .eq. icmv) ip2=0
152	if (iy1(ix ,jyp,indexh) .eq. icmv) ip3=0
153	if (iy1(ixp,jyp,indexh) .eq. icmv) ip4=0
154	ipsum= ip1+ip2+ip3+ip4
155	if (ipsum .eq. 0) then
156	yi1(m)=icmv
157	else
158	yi1(m)=(ip1p1iy1(ix ,jy ,indexh) &
159	+ ip2p2iy1(ixp,jy ,indexh) &
160	+ ip3p3iy1(ix ,jyp,indexh) &
161	+ ip4p4iy1(ixp,jyp,indexh))/ipsum
162	endif
163
164	ip1=1
165	ip2=1
166	ip3=1
167	ip4=1
168	if (iy2(ix ,jy ,indexh) .eq. icmv) ip1=0
169	if (iy2(ixp,jy ,indexh) .eq. icmv) ip2=0
170	if (iy2(ix ,jyp,indexh) .eq. icmv) ip3=0
171	if (iy2(ixp,jyp,indexh) .eq. icmv) ip4=0
172	ipsum= ip1+ip2+ip3+ip4
173	if (ipsum .eq. 0) then
174	yi2(m)=icmv
175	else
176	yi2(m)=(ip1p1iy2(ix ,jy ,indexh) &
177	+ ip2p2iy2(ixp,jy ,indexh) &
178	+ ip3p3iy2(ix ,jyp,indexh) &
179	+ ip4p4iy2(ixp,jyp,indexh))/ipsum
180	endif
181	!CPS end clouds
182
183	end do
184
185
186	!************************************
187	! 2.) Temporal interpolation (linear)
188	!************************************
189
190	if (abs(itime) .lt. abs(itime1)) then
191	print*,'interpol_rain.f90'
192	print*,itime,itime1,itime2
193	stop 'ITIME PROBLEM'
194	endif
195
196
197	dt1=real(itime-itime1)
198	dt2=real(itime2-itime)
199	dt=dt1+dt2
200
201	yint1=(y1(1)dt2+y1(2)dt1)/dt
202	yint2=(y2(1)dt2+y2(2)dt1)/dt
203	yint3=(y3(1)dt2+y3(2)dt1)/dt
204
205
206	!PS clouds:
207	intiy1=(yi1(1)dt2 + yi1(2)dt1)/dt
208	if (yi1(1) .eq. float(icmv)) intiy1=yi1(2)
209	if (yi1(2) .eq. float(icmv)) intiy1=yi1(1)
210
211	intiy2=(yi2(1)dt2 + yi2(2)dt1)/dt
212	if (yi2(1) .eq. float(icmv)) intiy2=yi2(2)
213	if (yi2(2) .eq. float(icmv)) intiy2=yi2(1)
214
215	if (intiy1 .ne. icmv .and. intiy2 .ne. icmv) then
216	intiy2 = intiy1 + intiy2 ! convert cloud thickness to cloud top
217	else
218	intiy1=icmv
219	intiy2=icmv
220	endif
221	!PS end clouds
222
223	end subroutine interpol_rain

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

Download in other formats:

Original Format