Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

interpol_rain.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 10 years ago
sources for flexwrf v3.1
File size: 10.0 KB

Line
1	!***********************************************************************
2	!* Copyright 2012,2013 *
3	!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine, *
4	!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5	!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso, *
6	!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann, *
7	!* *
8	!* This file is part of FLEXPART WRF *
9	!* *
10	!* FLEXPART is free software: you can redistribute it and/or modify *
11	!* it under the terms of the GNU General Public License as published by*
12	!* the Free Software Foundation, either version 3 of the License, or *
13	!* (at your option) any later version. *
14	!* *
15	!* FLEXPART is distributed in the hope that it will be useful, *
16	!* but WITHOUT ANY WARRANTY; without even the implied warranty of *
17	!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18	!* GNU General Public License for more details. *
19	!* *
20	!* You should have received a copy of the GNU General Public License *
21	!* along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
22	!***********************************************************************
23	subroutine interpol_rain(yy1,yy2,yy3,iy1,iy2,nxmax,nymax,nzmax,nx, &
24	ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, &
25	intiy1,intiy2,icmv)
26	! i i i i i i i
27	! i i i i i i i i o o o
28	!****************************************************************************
29	! *
30	! Interpolation of meteorological fields on 2-d model layers. *
31	! In horizontal direction bilinear interpolation interpolation is used. *
32	! Temporally a linear interpolation is used. *
33	! Three fields are interpolated at the same time. *
34	! *
35	! This is a special version of levlininterpol to save CPU time. *
36	! *
37	! 1 first time *
38	! 2 second time *
39	! *
40	! *
41	! Author: A. Stohl *
42	! *
43	! 30 August 1996 *
44	! *
45	!****************************************************************************
46	! *
47	! Variables: *
48	! *
49	! dt1,dt2 time differences between fields and current position *
50	! dz1,dz2 z distance between levels and current position *
51	! height(nzmax) heights of the model levels *
52	! indexh help variable *
53	! indz the level closest to the current trajectory position *
54	! indzh help variable *
55	! itime current time *
56	! itime1 time of the first wind field *
57	! itime2 time of the second wind field *
58	! ix,jy x,y coordinates of lower left subgrid point *
59	! level level at which interpolation shall be done *
60	! memind(3) points to the places of the wind fields *
61	! nx,ny actual field dimensions in x,y and z direction *
62	! nxmax,nymax,nzmax maximum field dimensions in x,y and z direction *
63	! xt current x coordinate *
64	! yint the final interpolated value *
65	! yt current y coordinate *
66	! yy(0:nxmax,0:nymax,nzmax,3) meteorological field used for interpolation *
67	! zt current z coordinate *
68	! *
69	! 19-Oct-2007 WANG, W change cumulated rain to hourly rain
70	! yy1 is accumulated largescale precipitation (mm)
71	! yy2 is ''''''''''''convective '''''''''''''''''
72	!****************************************************************************
73
74	implicit none
75
76	integer :: nx,ny,nxmax,nymax,nzmax,memind(2),m,ix,jy,ixp,jyp
77	! integer :: itime,itime1,itime2,level,indexh,i1,i2
78	integer :: itime,itime1,itime2,level,indexh,i1,i2,ip1,ip2,ip3,ip4
79	integer :: intiy1,intiy2,ipsum,icmv
80
81	real :: yy1(0:nxmax-1,0:nymax-1,nzmax,2)
82	real :: yy2(0:nxmax-1,0:nymax-1,nzmax,2)
83	real :: yy3(0:nxmax-1,0:nymax-1,nzmax,2)
84
85	integer :: iy1(0:nxmax-1,0:nymax-1,2),iy2(0:nxmax-1,0:nymax-1,2)
86
87	! real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
88	real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),yi1(2),yi2(2)
89
90	! real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
91	real :: xt,yt,yint1,yint2,yint3,yint4,p1,p2,p3,p4
92
93
94
95	! If point at border of grid -> small displacement into grid
96	!***********************************************************
97
98	if (xt.ge.real(nx-1)) xt=real(nx-1)-0.00001
99	if (yt.ge.real(ny-1)) yt=real(ny-1)-0.00001
100
101
102
103	!**********************************************************************
104	! 1.) Bilinear horizontal interpolation
105	! This has to be done separately for 2 fields (Temporal)
106	!*******************************************************
107
108	! Determine the lower left corner and its distance to the current position
109	!*************************************************************************
110
111	ix=int(xt)
112	jy=int(yt)
113	ixp=ix+1
114	jyp=jy+1
115	ddx=xt-real(ix)
116	ddy=yt-real(jy)
117	rddx=1.-ddx
118	rddy=1.-ddy
119	p1=rddx*rddy
120	p2=ddx*rddy
121	p3=rddx*ddy
122	p4=ddx*ddy
123
124
125	! Loop over 2 time steps
126	!***********************
127
128	! y1 and y2 are accumulated rain, need change to hourly rain
129
130	i1=memind(1)
131	i2=memind(2)
132	!! time interval between two fields, second to hour
133	dt=real(itime2-itime1)/3600.0
134
135
136	yint1=p1*(yy1(ix ,jy ,level,i2)- &
137	yy1(ix ,jy ,level,i1)) &
138	+ p2*(yy1(ixp,jy ,level,i2)- &
139	yy1(ixp,jy ,level,i1)) &
140	+ p3*(yy1(ix ,jyp,level,i2)- &
141	yy1(ix ,jyp,level,i1)) &
142	+ p4*(yy1(ixp,jyp,level,i2)- &
143	yy1(ixp,jyp,level,i1))
144	yint1=yint1/dt
145
146	yint2=p1*(yy2(ix ,jy ,level,i2)- &
147	yy2(ix ,jy ,level,i1)) &
148	+ p2*(yy2(ixp,jy ,level,i2)- &
149	yy2(ixp,jy ,level,i1)) &
150	+ p3*(yy2(ix ,jyp,level,i2)- &
151	yy2(ix ,jyp,level,i1)) &
152	+ p4*(yy2(ixp,jyp,level,i2)- &
153	yy2(ixp,jyp,level,i1))
154	yint2=yint2/dt
155
156
157	! Y3 is cloud fraction in an hour
158	do m=1,2
159	indexh=memind(m)
160	y3(m)=p1*yy3(ix ,jy ,level,indexh) &
161	+ p2*yy3(ixp,jy ,level,indexh) &
162	+ p3*yy3(ix ,jyp,level,indexh) &
163	+ p4*yy3(ixp,jyp,level,indexh)
164	enddo
165
166
167	! CDA new clouds
168
169	do m=1,2
170	indexh=memind(m)
171
172	ip1=1
173	ip2=1
174	ip3=1
175	ip4=1
176	if (iy1(ix ,jy ,indexh) .eq. icmv) ip1=0
177	if (iy1(ixp,jy ,indexh) .eq. icmv) ip2=0
178	if (iy1(ix ,jyp,indexh) .eq. icmv) ip3=0
179	if (iy1(ixp,jyp,indexh) .eq. icmv) ip4=0
180	ipsum= ip1+ip2+ip3+ip4
181	if (ipsum .eq. 0) then
182	yi1(m)=icmv
183	else
184	yi1(m)=(ip1p1iy1(ix ,jy ,indexh) &
185	+ ip2p2iy1(ixp,jy ,indexh) &
186	+ ip3p3iy1(ix ,jyp,indexh) &
187	+ ip4p4iy1(ixp,jyp,indexh))/ipsum
188	endif
189
190	ip1=1
191	ip2=1
192	ip3=1
193	ip4=1
194	if (iy2(ix ,jy ,indexh) .eq. icmv) ip1=0
195	if (iy2(ixp,jy ,indexh) .eq. icmv) ip2=0
196	if (iy2(ix ,jyp,indexh) .eq. icmv) ip3=0
197	if (iy2(ixp,jyp,indexh) .eq. icmv) ip4=0
198	ipsum= ip1+ip2+ip3+ip4
199	if (ipsum .eq. 0) then
200	yi2(m)=icmv
201	else
202	yi2(m)=(ip1p1iy2(ix ,jy ,indexh) &
203	+ ip2p2iy2(ixp,jy ,indexh) &
204	+ ip3p3iy2(ix ,jyp,indexh) &
205	+ ip4p4iy2(ixp,jyp,indexh))/ipsum
206	endif
207	enddo
208	!CPS end clouds
209
210
211	10 continue
212
213
214	!************************************
215	! 2.) Temporal interpolation (linear)
216	!************************************
217
218	dt1=real(itime-itime1)
219	dt2=real(itime2-itime)
220	dt=dt1+dt2
221
222	! yint1=(y1(1)dt2+y1(2)dt1)/dt
223	! yint2=(y2(1)dt2+y2(2)dt1)/dt
224	yint3=(y3(1)dt2+y3(2)dt1)/dt
225
226
227	!CPS clouds:
228	intiy1=(yi1(1)dt2 + yi1(2)dt1)/dt
229	if (yi1(1) .eq. float(icmv)) intiy1=yi1(2)
230	if (yi1(2) .eq. float(icmv)) intiy1=yi1(1)
231
232	intiy2=(yi2(1)dt2 + yi2(2)dt1)/dt
233	if (yi2(1) .eq. float(icmv)) intiy2=yi2(2)
234	if (yi2(2) .eq. float(icmv)) intiy2=yi2(1)
235
236	if (intiy1 .ne. icmv .and. intiy2 .ne. icmv) then
237	intiy2 = intiy1 + intiy2 ! convert cloud thickness to cloud top
238	else
239	intiy1=icmv
240	intiy2=icmv
241	endif
242	!CPS end clouds
243
244
245
246
247
248	end subroutine interpol_rain
249

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

Download in other formats:

Original Format