Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

interpol_rain_nests.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 11 years ago
sources for flexwrf v3.1
File size: 10.9 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_nests(yy1,yy2,yy3,iy1,iy2,nxmaxn,nymaxn,nzmax, &
24	maxnests,ngrid,nxn,nyn,memind,xt,yt,level,itime1,itime2,itime, &
25	yint1,yint2,yint3,intiy1,intiy2,icmv)
26	! i i i i i i
27	! i i i i i i i i i i i
28	! o o o
29	!****************************************************************************
30	! *
31	! Interpolation of meteorological fields on 2-d model layers for nested *
32	! grids. This routine is related to levlin3interpol.f for the mother domain*
33	! *
34	! In horizontal direction bilinear interpolation interpolation is used. *
35	! Temporally a linear interpolation is used. *
36	! Three fields are interpolated at the same time. *
37	! *
38	! This is a special version of levlininterpol to save CPU time. *
39	! *
40	! 1 first time *
41	! 2 second time *
42	! *
43	! *
44	! Author: A. Stohl *
45	! *
46	! 15 March 2000 *
47	! *
48	!****************************************************************************
49	! *
50	! Variables: *
51	! *
52	! dt1,dt2 time differences between fields and current position *
53	! dz1,dz2 z distance between levels and current position *
54	! height(nzmax) heights of the model levels *
55	! indexh help variable *
56	! indz the level closest to the current trajectory position *
57	! indzh help variable *
58	! itime current time *
59	! itime1 time of the first wind field *
60	! itime2 time of the second wind field *
61	! ix,jy x,y coordinates of lower left subgrid point *
62	! level level at which interpolation shall be done *
63	! memind(3) points to the places of the wind fields *
64	! nx,ny actual field dimensions in x,y and z direction *
65	! nxmax,nymax,nzmax maximum field dimensions in x,y and z direction *
66	! xt current x coordinate *
67	! yint the final interpolated value *
68	! yt current y coordinate *
69	! yy(0:nxmax,0:nymax,nzmax,3) meteorological field used for interpolation *
70	! zt current z coordinate *
71	! *
72	! Changed 10/22/2007 yy1,yy2 are accumulated rain (mm)
73	! convert them into hourly rain (mm/hr)
74	!****************************************************************************
75
76	implicit none
77
78	integer :: maxnests,ngrid
79	integer :: nxn(maxnests),nyn(maxnests),nxmaxn,nymaxn,nzmax,memind(2)
80	integer :: m,ix,jy,ixp,jyp,itime,itime1,itime2,level,indexh,i1,i2
81	integer :: ip1,ip2,ip3,ip4
82	integer :: intiy1,intiy2,ipsum,icmv
83	real :: yy1(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
84	real :: yy2(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
85	real :: yy3(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
86	integer :: iy1(0:nxmaxn-1,0:nymaxn-1,2,maxnests)
87	integer :: iy2(0:nxmaxn-1,0:nymaxn-1,2,maxnests)
88	real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),yi1(2),yi2(2)
89	real :: xt,yt,yint1,yint2,yint3,yint4,p1,p2,p3,p4
90
91
92
93	! If point at border of grid -> small displacement into grid
94	!***********************************************************
95
96	if (xt.ge.real(nxn(ngrid)-1)) xt=real(nxn(ngrid)-1)-0.00001
97	if (yt.ge.real(nyn(ngrid)-1)) yt=real(nyn(ngrid)-1)-0.00001
98
99
100
101	!**********************************************************************
102	! 1.) Bilinear horizontal interpolation
103	! This has to be done separately for 2 fields (Temporal)
104	!*******************************************************
105
106	! Determine the lower left corner and its distance to the current position
107	!*************************************************************************
108
109	ix=int(xt)
110	jy=int(yt)
111	ixp=ix+1
112	jyp=jy+1
113	ddx=xt-real(ix)
114	ddy=yt-real(jy)
115	rddx=1.-ddx
116	rddy=1.-ddy
117	p1=rddx*rddy
118	p2=ddx*rddy
119	p3=rddx*ddy
120	p4=ddx*ddy
121
122
123	! Loop over 2 time steps
124	!***********************
125
126	! y1 and y2 are accumulated rain, need change to hourly rain
127
128	i1=memind(1)
129	i2=memind(2)
130	! time interval between two fields, second to hour
131	dt=real(itime2-itime1)/3600.0
132
133	yint1=p1*(yy1(ix ,jy ,level,i2,ngrid)- &
134	yy1(ix ,jy ,level,i1,ngrid)) &
135	+ p2*(yy1(ixp,jy ,level,i2,ngrid)- &
136	yy1(ixp,jy ,level,i1,ngrid)) &
137	+ p3*(yy1(ix ,jyp,level,i2,ngrid)- &
138	yy1(ix ,jyp,level,i1,ngrid)) &
139	+ p4*(yy1(ixp,jyp,level,i2,ngrid)- &
140	yy1(ixp,jyp,level,i1,ngrid))
141	yint1=yint1/dt
142
143	yint2=p1*(yy2(ix ,jy ,level,i2,ngrid)- &
144	yy2(ix ,jy ,level,i1,ngrid)) &
145	+ p2*(yy2(ixp,jy ,level,i2,ngrid)- &
146	yy2(ixp,jy ,level,i1,ngrid)) &
147	+ p3*(yy2(ix ,jyp,level,i2,ngrid)- &
148	yy2(ix ,jyp,level,i1,ngrid)) &
149	+ p4*(yy2(ixp,jyp,level,i2,ngrid)- &
150	yy2(ixp,jyp,level,i1,ngrid))
151	yint2=yint2/dt
152
153
154	! Y3 is cloud fraction in an hour
155
156
157	do m=1,2
158	indexh=memind(m)
159
160	y3(m)=p1*yy3(ix ,jy ,level,indexh,ngrid) &
161	+ p2*yy3(ixp,jy ,level,indexh,ngrid) &
162	+ p3*yy3(ix ,jyp,level,indexh,ngrid) &
163	+ p4*yy3(ixp,jyp,level,indexh,ngrid)
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,ngrid) .eq. icmv) ip1=0
177	if (iy1(ixp,jy ,indexh,ngrid) .eq. icmv) ip2=0
178	if (iy1(ix ,jyp,indexh,ngrid) .eq. icmv) ip3=0
179	if (iy1(ixp,jyp,indexh,ngrid) .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,ngrid) &
185	+ ip2p2iy1(ixp,jy ,indexh,ngrid) &
186	+ ip3p3iy1(ix ,jyp,indexh,ngrid) &
187	+ ip4p4iy1(ixp,jyp,indexh,ngrid))/ipsum
188	endif
189
190	ip1=1
191	ip2=1
192	ip3=1
193	ip4=1
194	if (iy2(ix ,jy ,indexh,ngrid) .eq. icmv) ip1=0
195	if (iy2(ixp,jy ,indexh,ngrid) .eq. icmv) ip2=0
196	if (iy2(ix ,jyp,indexh,ngrid) .eq. icmv) ip3=0
197	if (iy2(ixp,jyp,indexh,ngrid) .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,ngrid) &
203	+ ip2p2iy2(ixp,jy ,indexh,ngrid) &
204	+ ip3p3iy2(ix ,jyp,indexh,ngrid) &
205	+ ip4p4iy2(ixp,jyp,indexh,ngrid))/ipsum
206	endif
207	ip1=1
208	ip2=1
209	ip3=1
210	ip4=1
211	if (iy2(ix ,jy ,indexh,ngrid) .eq. icmv) ip1=0
212	if (iy2(ixp,jy ,indexh,ngrid) .eq. icmv) ip2=0
213	if (iy2(ix ,jyp,indexh,ngrid) .eq. icmv) ip3=0
214	if (iy2(ixp,jyp,indexh,ngrid) .eq. icmv) ip4=0
215	ipsum= ip1+ip2+ip3+ip4
216	if (ipsum .eq. 0) then
217	yi2(m)=icmv
218	else
219	yi2(m)=(ip1p1iy2(ix ,jy ,indexh,ngrid) &
220	+ ip2p2iy2(ixp,jy ,indexh,ngrid) &
221	+ ip3p3iy2(ix ,jyp,indexh,ngrid) &
222	+ ip4p4iy2(ixp,jyp,indexh,ngrid))/ipsum
223	endif
224	enddo
225	!CPS end clouds
226
227
228	10 continue
229
230
231	!************************************
232	! 2.) Temporal interpolation (linear)
233	!************************************
234
235	dt1=real(itime-itime1)
236	dt2=real(itime2-itime)
237	dt=dt1+dt2
238
239	! yint1=(y1(1)dt2+y1(2)dt1)/dt
240	! yint2=(y2(1)dt2+y2(2)dt1)/dt
241	yint3=(y3(1)dt2+y3(2)dt1)/dt
242
243
244	!CPS clouds:
245	intiy1=(yi1(1)dt2 + yi1(2)dt1)/dt
246	if (yi1(1) .eq. float(icmv)) intiy1=yi1(2)
247	if (yi1(2) .eq. float(icmv)) intiy1=yi1(1)
248
249	intiy2=(yi2(1)dt2 + yi2(2)dt1)/dt
250	if (yi2(1) .eq. float(icmv)) intiy2=yi2(2)
251	if (yi2(2) .eq. float(icmv)) intiy2=yi2(1)
252
253	if (intiy1 .ne. icmv .and. intiy2 .ne. icmv) then
254	intiy2 = intiy1 + intiy2 ! convert cloud thickness to cloud top
255	else
256	intiy1=icmv
257	intiy2=icmv
258	endif
259	!CPS end clouds
260
261
262	end subroutine interpol_rain_nests
263

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

Download in other formats:

Original Format