Context Navigation

interpol_rain_nests.f90 @ 37

Last change on this file since 37 was 37, checked in by pesei, 9 years ago
Wet dep quick fix and other small changes. Wet depo quick fix not final yet.
File size: 9.5 KB

Line
1	!**********************************************************************
2	! Copyright 1998-2015 *
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_nests(yy1,yy2,yy3,nxmaxn,nymaxn,nzmax, &
23	maxnests,ngrid,nxn,nyn, &
24	memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, &
25	icbot,ictop)
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	! Petra Seibert, 2011/2012-2015:
50	! Add interpolation of cloud bottom and cloud thickness
51	! for fix to SE's new wet scavenging scheme
52	! *
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	use com_mod, only: icloudbotn, icloudthckn
80	use par_mod, only: icmv
81
82	implicit none
83
84	integer :: maxnests,ngrid
85	integer :: nxn(maxnests),nyn(maxnests),nxmaxn,nymaxn,nzmax,memind(2),m
86	integer :: ix,jy,ixp,jyp
87	integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4
88	integer :: icbot,ictop,icthck,ipsum
89	real :: yy1(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
90	real :: yy2(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
91	real :: yy3(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
92	real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),cbot(2),cthck(2)
93	real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
94
95	! If point at border of grid -> small displacement into grid
96	!***********************************************************
97
98	if (xt.ge.(real(nxn(ngrid)-1)-0.0001)) xt=real(nxn(ngrid)-1)-0.0001
99	if (yt.ge.(real(nyn(ngrid)-1)-0.0001)) yt=real(nyn(ngrid)-1)-0.0001
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	memind_loop: &
129	do m=1,2
130
131	indexh=memind(m)
132
133	y1(m)=p1*yy1(ix ,jy ,level,indexh,ngrid) &
134	+ p2*yy1(ixp,jy ,level,indexh,ngrid) &
135	+ p3*yy1(ix ,jyp,level,indexh,ngrid) &
136	+ p4*yy1(ixp,jyp,level,indexh,ngrid)
137	y2(m)=p1*yy2(ix ,jy ,level,indexh,ngrid) &
138	+ p2*yy2(ixp,jy ,level,indexh,ngrid) &
139	+ p3*yy2(ix ,jyp,level,indexh,ngrid) &
140	+ p4*yy2(ixp,jyp,level,indexh,ngrid)
141	y3(m)=p1*yy3(ix ,jy ,level,indexh,ngrid) &
142	+ p2*yy3(ixp,jy ,level,indexh,ngrid) &
143	+ p3*yy3(ix ,jyp,level,indexh,ngrid) &
144	+ p4*yy3(ixp,jyp,level,indexh,ngrid)
145
146	! PS clouds:
147	ip1=1
148	ip2=1
149	ip3=1
150	ip4=1
151	if (icloudbotn(ix ,jy ,indexh,ngrid) .eq. icmv) ip1=0
152	if (icloudbotn(ixp,jy ,indexh,ngrid) .eq. icmv) ip2=0
153	if (icloudbotn(ix ,jyp,indexh,ngrid) .eq. icmv) ip3=0
154	if (icloudbotn(ixp,jyp,indexh,ngrid) .eq. icmv) ip4=0
155	ipsum = ip1+ip2+ip3+ip4
156	if (ipsum .eq. 0) then
157	cbot(m)=icmv
158	else
159	cbot(m)=(ip1p1icloudbotn(ix ,jy ,indexh,ngrid) &
160	+ ip2p2icloudbotn(ixp,jy ,indexh,ngrid) &
161	+ ip3p3icloudbotn(ix ,jyp,indexh,ngrid) &
162	+ ip4p4icloudbotn(ixp,jyp,indexh,ngrid)) / ipsum
163	endif
164
165	ip1=1
166	ip2=1
167	ip3=1
168	ip4=1
169	if (icloudthckn(ix ,jy ,indexh,ngrid) .eq. icmv) ip1=0
170	if (icloudthckn(ixp,jy ,indexh,ngrid) .eq. icmv) ip2=0
171	if (icloudthckn(ix ,jyp,indexh,ngrid) .eq. icmv) ip3=0
172	if (icloudthckn(ixp,jyp,indexh,ngrid) .eq. icmv) ip4=0
173	ipsum = ip1+ip2+ip3+ip4
174	if (ipsum .eq. 0) then
175	cthck(m)=icmv
176	else
177	cthck(m)=(ip1p1icloudthckn(ix ,jy ,indexh,ngrid) &
178	+ ip2p2icloudthckn(ixp,jy ,indexh,ngrid) &
179	+ ip3p3icloudthckn(ix ,jyp,indexh,ngrid) &
180	+ ip4p4icloudthckn(ixp,jyp,indexh,ngrid)) / ipsum
181	endif
182	! PS end clouds
183
184	enddo memind_loop
185
186
187	!************************************
188	! 2.) Temporal interpolation (linear)
189	!************************************
190
191	dt1=real(itime-itime1)
192	dt2=real(itime2-itime)
193	dt=dt1+dt2
194
195	yint1=(y1(1)dt2+y1(2)dt1)/dt
196	yint2=(y2(1)dt2+y2(2)dt1)/dt
197	yint3=(y3(1)dt2+y3(2)dt1)/dt
198
199
200	! PS clouds:
201	icbot = nint( (cbot(1)dt2 + cbot(2)dt1)/dt )
202	if (nint(cbot(1)) .eq. icmv) icbot=cbot(2)
203	if (nint(cbot(2)) .eq. icmv) icbot=cbot(1)
204
205	icthck = nint( (cthck(1)dt2 + cthck(2)dt1)/dt )
206	if (nint(cthck(1)) .eq. icmv) icthck=cthck(2)
207	if (nint(cthck(2)) .eq. icmv) icthck=cthck(1)
208
209	if (icbot .ne. icmv .and. icthck .ne. icmv) then
210	ictop = icbot + icthck ! convert cloud thickness to cloud top
211	else
212	icbot=icmv
213	ictop=icmv
214	endif
215	! PS end clouds
216
217
218	end subroutine interpol_rain_nests

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

Download in other formats:

Original Format