source: trunk/src/interpol_rain.f90 @ 20

Last change on this file since 20 was 20, checked in by igpis, 9 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
RevLine 
[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)=(ip1*p1*iy1(ix ,jy ,indexh)    &
159             + ip2*p2*iy1(ixp,jy ,indexh)       &
160             + ip3*p3*iy1(ix ,jyp,indexh)       &
161             + ip4*p4*iy1(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)=(ip1*p1*iy2(ix ,jy ,indexh)  &
177             + ip2*p2*iy2(ixp,jy ,indexh)     &
178             + ip3*p3*iy2(ix ,jyp,indexh)     &
179             + ip4*p4*iy2(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]223end subroutine interpol_rain
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG