Context Navigation

interpol_rain.f90 @ 20

Last change on this file since 20 was 20, checked in by igpis, 10 years ago
move version 9.1.8 form branches to trunk. Contributions from HSO, saeck, pesei, NIK, RT, XKF, IP and others
File size: 9.4 KB

Rev	Line
[4]	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
[20]	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
[4]	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 *
[20]	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 *
[4]	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
[20]	82	!integer :: itime,itime1,itime2,level,indexh
	83	integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4
	84	integer :: intiy1,intiy2,ipsum,icmv
[4]	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)
[20]	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
[4]	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)
[20]	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
[4]	183	end do
	184
	185
	186	!************************************
	187	! 2.) Temporal interpolation (linear)
	188	!************************************
	189
[20]	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
[4]	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
[20]	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
[4]	223	end subroutine interpol_rain

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

Download in other formats:

Original Format