interpol_rain_nests.f90 6.39 KB
Newer Older
1
2
! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
! SPDX-License-Identifier: GPL-3.0-or-later
3

Matthias Langer's avatar
 
Matthias Langer committed
4
subroutine interpol_rain_nests(yy1,yy2,yy3,nxmaxn,nymaxn,nzmax, &
5
       maxnests,ngrid,nxn,nyn,iwftouse,xt,yt,level,itime1,itime2,itime, &
Matthias Langer's avatar
 
Matthias Langer committed
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
       yint1,yint2,yint3)
  !                                i   i   i    i      i      i
  !   i       i    i   i    i    i  i    i     i      i      i
  !  o     o     o
  !****************************************************************************
  !                                                                           *
  !  Interpolation of meteorological fields on 2-d model layers for nested    *
  !  grids. This routine is related to levlin3interpol.f for the mother domain*
  !                                                                           *
  !  In horizontal direction bilinear interpolation interpolation is used.    *
  !  Temporally a linear interpolation is used.                               *
  !  Three fields are interpolated at the same time.                          *
  !                                                                           *
  !  This is a special version of levlininterpol to save CPU time.            *
  !                                                                           *
  !  1 first time                                                             *
  !  2 second time                                                            *
  !                                                                           *
  !                                                                           *
  !     Author: A. Stohl                                                      *
  !                                                                           *
  !     15 March 2000                                                         *
  !                                                                           *
  !****************************************************************************
  !                                                                           *
  ! Variables:                                                                *
  !                                                                           *
  ! dt1,dt2              time differences between fields and current position *
  ! dz1,dz2              z distance between levels and current position       *
  ! height(nzmax)        heights of the model levels                          *
  ! indexh               help variable                                        *
  ! indz                 the level closest to the current trajectory position *
  ! indzh                help variable                                        *
  ! itime                current time                                         *
  ! itime1               time of the first wind field                         *
  ! itime2               time of the second wind field                        *
  ! ix,jy                x,y coordinates of lower left subgrid point          *
  ! level                level at which interpolation shall be done           *
44
  ! iwftouse             points to the place of the wind field                *
Matthias Langer's avatar
 
Matthias Langer committed
45
46
47
48
49
50
51
52
53
  ! nx,ny                actual field dimensions in x,y and z direction       *
  ! nxmax,nymax,nzmax    maximum field dimensions in x,y and z direction      *
  ! xt                   current x coordinate                                 *
  ! yint                 the final interpolated value                         *
  ! yt                   current y coordinate                                 *
  ! yy(0:nxmax,0:nymax,nzmax,3) meteorological field used for interpolation   *
  ! zt                   current z coordinate                                 *
  !                                                                           *
  !****************************************************************************
54
  use par_mod, only: numwfmem
Matthias Langer's avatar
 
Matthias Langer committed
55
56
57
58

  implicit none

  integer :: maxnests,ngrid
59
  integer :: nxn(maxnests),nyn(maxnests),nxmaxn,nymaxn,nzmax,iwftouse
Matthias Langer's avatar
 
Matthias Langer committed
60
  integer :: m,ix,jy,ixp,jyp,itime,itime1,itime2,level,indexh
61
62
63
  real :: yy1(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,maxnests)
  real :: yy2(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,maxnests)
  real :: yy3(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,maxnests)
Matthias Langer's avatar
 
Matthias Langer committed
64
65
66
67
68
69
70
71
  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
  real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4



  ! If point at border of grid -> small displacement into grid
  !***********************************************************

72
73
74
75
76
77
78
79
  ! if (xt.ge.(real(nxn(ngrid)-1)-0.0001)) &
  !      xt=real(nxn(ngrid)-1)-0.0001
  ! if (yt.ge.(real(nyn(ngrid)-1)-0.0001)) &
  !      yt=real(nyn(ngrid)-1)-0.0001

! ESO make it consistent with interpol_rain
  if (xt.ge.(real(nxn(ngrid)-1))) xt=real(nxn(ngrid)-1)-0.00001
  if (yt.ge.(real(nyn(ngrid)-1))) yt=real(nyn(ngrid)-1)-0.00001
Matthias Langer's avatar
 
Matthias Langer committed
80
81
82
83
84
85
86
87
88
89
90
91
92



  !**********************************************************************
  ! 1.) Bilinear horizontal interpolation
  ! This has to be done separately for 2 fields (Temporal)
  !*******************************************************

  ! Determine the lower left corner and its distance to the current position
  !*************************************************************************

  ix=int(xt)
  jy=int(yt)
93

Matthias Langer's avatar
 
Matthias Langer committed
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
  ixp=ix+1
  jyp=jy+1
  ddx=xt-real(ix)
  ddy=yt-real(jy)
  rddx=1.-ddx
  rddy=1.-ddy
  p1=rddx*rddy
  p2=ddx*rddy
  p3=rddx*ddy
  p4=ddx*ddy


  ! Loop over 2 time steps
  !***********************

109
110
111
!  do m=1,2
!    indexh=memind(m)
    indexh=iwftouse
Matthias Langer's avatar
 
Matthias Langer committed
112

113
    y1(1)=p1*yy1(ix ,jy ,level,indexh,ngrid) &
Matthias Langer's avatar
 
Matthias Langer committed
114
115
116
         + p2*yy1(ixp,jy ,level,indexh,ngrid) &
         + p3*yy1(ix ,jyp,level,indexh,ngrid) &
         + p4*yy1(ixp,jyp,level,indexh,ngrid)
117
    y2(1)=p1*yy2(ix ,jy ,level,indexh,ngrid) &
Matthias Langer's avatar
 
Matthias Langer committed
118
119
120
         + p2*yy2(ixp,jy ,level,indexh,ngrid) &
         + p3*yy2(ix ,jyp,level,indexh,ngrid) &
         + p4*yy2(ixp,jyp,level,indexh,ngrid)
121
    y3(1)=p1*yy3(ix ,jy ,level,indexh,ngrid) &
Matthias Langer's avatar
 
Matthias Langer committed
122
123
124
         + p2*yy3(ixp,jy ,level,indexh,ngrid) &
         + p3*yy3(ix ,jyp,level,indexh,ngrid) &
         + p4*yy3(ixp,jyp,level,indexh,ngrid)
125
!  end do
Matthias Langer's avatar
 
Matthias Langer committed
126
127
128
129
130
131


  !************************************
  ! 2.) Temporal interpolation (linear)
  !************************************

132
133
134
  ! dt1=real(itime-itime1)
  ! dt2=real(itime2-itime)
  ! dt=dt1+dt2
Matthias Langer's avatar
 
Matthias Langer committed
135

136
137
138
  ! yint1=(y1(1)*dt2+y1(2)*dt1)/dt
  ! yint2=(y2(1)*dt2+y2(2)*dt1)/dt
  ! yint3=(y3(1)*dt2+y3(2)*dt1)/dt
Matthias Langer's avatar
 
Matthias Langer committed
139

140
141
142
   yint1=y1(1)
   yint2=y2(1)
   yint3=y3(1)
Matthias Langer's avatar
 
Matthias Langer committed
143
144

end subroutine interpol_rain_nests