1
2
3
4
5
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
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
|
subroutine intermod(xt,yt,zt,indexf,itime,init,firststep,
+iter,lkindz,ngrid,uint,vint)
C i i i i i i i
C i i i o o
********************************************************************************
* *
* Interpolation routine for trajectories on model layers. *
* *
* Author: A. Stohl *
* *
* 27 April 1994 *
* *
* Update: 23 December 1998 (Use of global domain and nesting) *
* *
********************************************************************************
* *
* Variables: *
* firstep .true. for first iteration of petterssen *
* idiff [s] Temporal distance between the windfields used for interpol*
* init .true. for first time step of trajectory *
* iter number of iteration step *
* itime [s] current temporal position *
* lkindz unit of z coordinate (1:masl, 2:magl, 3:hPa) *
* memtime(3) [s] times of the wind fields in memory *
* ngrid index which grid is to be used *
* nstop =greater 0, if trajectory calculation is finished *
* uint,vint [m/s] interpolated wind components *
* xt,yt,zt coordinates position for which wind data shall be calculat*
* *
********************************************************************************
include 'includepar'
include 'includecom'
integer itime,indexf,iter,lkindz,ngrid
real xt,yt,zt,uint,vint,psint,eta,xtn,ytn
logical firststep,init
C Determine nested grid coordinates
***********************************
if (ngrid.gt.0) then
xtn=(xt-xln(ngrid))*xresoln(ngrid)
ytn=(yt-yln(ngrid))*yresoln(ngrid)
endif
C For a new trajectory and the first iteration of petterssen:
C Transformation from height in [m] coordinate to height in eta coordinate.
C Or: Transformation from hPa to eta.
***************************************************************************
if (init.and.firststep) then
C Calculate the surface pressure.
*********************************
if (ngrid.gt.0) then ! nested grid
call levlininterpoln(psn,maxnests,nxmaxn,nymaxn,1,ngrid,
+ nxn,nyn,memind,xtn,ytn,1,memtime(indexf),memtime(indexf+1),
+ itime,indexf,psint)
else ! mother grid
call levlininterpol(ps,nxmax,nymax,1,nx,ny,memind,xt,yt,
+ 1,memtime(indexf),memtime(indexf+1),itime,indexf,psint)
endif
if ((lkindz.eq.1).or.(lkindz.eq.2)) then
if (ngrid.gt.0) then ! nested grid
call etatrafo(xtn,ytn,zt,memtime(indexf),memtime(indexf+1),
+ itime,indexf,ngrid,psint)
else ! mother grid
call etatrafo(xt,yt,zt,memtime(indexf),memtime(indexf+1),
+ itime,indexf,ngrid,psint)
endif
else if (lkindz.eq.3) then
zt=eta(psint,zt)
endif
endif
C Interpolation of wind field data is done.
C Either linear interpolation (if selected and for first iteration
C of petterssen) or bicubic interpolation
****************************************************************************
if (ngrid.eq.0) then ! mother grid
if ((inpolkind.eq.1).and.(iter.ne.1)) then
call interpol(uu,nxmax,nymax,nuvzmax,nx,ny,nuvz,memind,uvheight,
+ xt,yt,zt,memtime(indexf),memtime(indexf+1),itime,indexf,uint)
call interpol(vv,nxmax,nymax,nuvzmax,nx,ny,nuvz,memind,uvheight,
+ xt,yt,zt,memtime(indexf),memtime(indexf+1),itime,indexf,vint)
else ! bilinear interpolation
call lininterpol(uu,nxmax,nymax,nuvzmax,nx,ny,nuvz,memind,
+ uvheight,xt,yt,zt,memtime(indexf),memtime(indexf+1),itime,
+ indexf,uint)
call lininterpol(vv,nxmax,nymax,nuvzmax,nx,ny,nuvz,memind,
+ uvheight,xt,yt,zt,memtime(indexf),memtime(indexf+1),itime,
+ indexf,vint)
endif
else if (ngrid.gt.0) then ! nested grid
if ((inpolkind.eq.1).and.(iter.ne.1)) then
call interpoln(uun,maxnests,nxmaxn,nymaxn,nuvzmax,
+ ngrid,nxn,nyn,nuvz,memind,uvheight,xtn,ytn,zt,memtime(indexf),
+ memtime(indexf+1),itime,indexf,uint)
call interpoln(vvn,maxnests,nxmaxn,nymaxn,nuvzmax,
+ ngrid,nxn,nyn,nuvz,memind,uvheight,xtn,ytn,zt,memtime(indexf),
+ memtime(indexf+1),itime,indexf,vint)
else ! bilinear interpolation
call lininterpoln(uun,maxnests,nxmaxn,nymaxn,nuvzmax,
+ ngrid,nxn,nyn,nuvz,memind,uvheight,xtn,ytn,zt,memtime(indexf),
+ memtime(indexf+1),itime,indexf,uint)
call lininterpoln(vvn,maxnests,nxmaxn,nymaxn,nuvzmax,
+ ngrid,nxn,nyn,nuvz,memind,uvheight,xtn,ytn,zt,memtime(indexf),
+ memtime(indexf+1),itime,indexf,vint)
endif
else ! polar stereographic grid
! only linear interpolation used
call lininterpol(uupol,nxmax,nymax,nuvzmax,nx,ny,nuvz,memind,
+ uvheight,xt,yt,zt,memtime(indexf),memtime(indexf+1),itime,
+ indexf,uint)
call lininterpol(vvpol,nxmax,nymax,nuvzmax,nx,ny,nuvz,memind,
+ uvheight,xt,yt,zt,memtime(indexf),memtime(indexf+1),itime,
+ indexf,vint)
endif
return
end
|
