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
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
|
subroutine interisentrop(xt,yt,zt,thetaconst,indexf,itime,init,
+firststep,iter,lkindz,ngrid,uint,vint)
C i i i/o i/o i i i
C i i i i o o
********************************************************************************
* *
* Interpolation routine for isentropic trajectories. *
* *
* The user gives the vertical coordinate in metres. In the first time step of *
* the trajectory, metres are transformed to Kelvin. In all time steps these *
* Kelvin are kept constant and transformed onward to the eta coordinate. *
* *
* *
* Author: A. Stohl *
* *
* 27 April 1994 *
* *
* Update: 23 December 1998 (Use of global domain and nesting) *
* *
********************************************************************************
* *
* Variables: *
* firststep .true. for first step of petterssen, .false. for iteration*
* idiff [s] Temporal distance between the windfields used for interpol*
* indz1,indz2 indices of boundary of unstable region *
* init .true. for first time step of trajectory, .false. for othe*
* 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 *
* psint [Pa] interpolated surface pressure *
* temp [K] temperature at trajectory position *
* thetaconst [Pa] theta level, for which trajectories shall be calculated *
* uint,vint [m/s] interpolated wind components *
* uvheight heights in which u,v and t are given *
* xt,yt,zt coordinates position for which wind data shall be calculat*
* *
* Constants: *
* kappa exponent for calculating potential temperature *
* *
* Function: *
* pp calculates the pressure at a given eta level *
* *
********************************************************************************
include 'includepar'
include 'includecom'
integer itime,indexf,iter,indz1,indz2,i,lkindz,ngrid
real xt,yt,zt,uint,vint,psint,thetaconst,phelp,pp,xtn,ytn
real weight,uhelp,vhelp,qhelp,ph1,ph2,eta
logical init,firststep
indz1=0
indz2=0
C Determine nested grid coordinates
***********************************
if (ngrid.gt.0) then
xtn=(xt-xln(ngrid))*xresoln(ngrid)
ytn=(yt-yln(ngrid))*yresoln(ngrid)
endif
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
******************************************************************************
C For a new trajectory and the first iteration of petterssen:
C thetaconst is given in [m]
C 1. Transformation from height in [m] coordinate to height in eta coordinate.
C Or: Transformation from hPa to eta.
C 2. Transformation from height in eta to height in theta [K] coordinate.
******************************************************************************
if (init.and.firststep) then
if ((lkindz.eq.1).or.(lkindz.eq.2)) then
if (ngrid.gt.0) then ! nested grid
call etatrafo(xtn,ytn,thetaconst,memtime(indexf),
+ memtime(indexf+1),itime,indexf,ngrid,psint)
else ! mother grid
call etatrafo(xt,yt,thetaconst,memtime(indexf),
+ memtime(indexf+1),itime,indexf,ngrid,psint)
endif
zt=thetaconst ! now thetaconst is given in eta
else if (lkindz.eq.3) then
zt=eta(psint,zt)
endif
C Calculate pressure and potential temperature for current position
*******************************************************************
phelp=pp(psint,zt)
if (ngrid.gt.0) then ! nested grid
call lininterpoln(thn,maxnests,nxmaxn,nymaxn,nuvzmax,
+ ngrid,nxn,nyn,nuvz,memind,uvheight,xtn,ytn,zt,memtime(indexf),
+ memtime(indexf+1),itime,indexf,thetaconst)
else ! mother grid
call lininterpol(th,nxmax,nymax,nuvzmax,nx,ny,nuvz,memind,
+ uvheight,xt,yt,zt,memtime(indexf),memtime(indexf+1),itime,
+ indexf,thetaconst)
endif
else
********************************************************************************
C For all other time steps:
C Transformation from height in theta [K] coordinate to height in eta coordinate
********************************************************************************
if (ngrid.gt.0) then ! nested grid
call geteta(xtn,ytn,zt,thetaconst,memtime(indexf),
+ memtime(indexf+1),itime,indexf,psint,ngrid,indz1,indz2)
else ! mother grid
call geteta(xt,yt,zt,thetaconst,memtime(indexf),
+ memtime(indexf+1),itime,indexf,psint,ngrid,indz1,indz2)
endif
endif
C Interpolation of wind field data is done.
*****************************************************************************
if (indz1.eq.indz2) then ! stable level
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
else ! unstable layer
C -> take the layer averaged wind velocity (weighted with pressure increments)
******************************************************************************
weight=0.
uhelp=0.
vhelp=0.
qhelp=0.
ph1=akz(indz1)+bkz(indz1)*psint
do 10 i=indz1,indz2-1
if (ngrid.ge.0) then ! mother grid
if ((inpolkind.eq.1).and.(iter.ne.1)) then
call levinterpol(uu,nxmax,nymax,nuvzmax,nx,ny,memind,
+ xt,yt,i,memtime(indexf),memtime(indexf+1),itime,indexf,
+ uint)
call levinterpol(vv,nxmax,nymax,nuvzmax,nx,ny,memind,
+ xt,yt,i,memtime(indexf),memtime(indexf+1),itime,indexf,
+ vint)
else ! bilinear interpolation
call levlininterpol(uu,nxmax,nymax,nuvzmax,nx,ny,memind,
+ xt,yt,i,memtime(indexf),memtime(indexf+1),itime,indexf,
+ uint)
call levlininterpol(vv,nxmax,nymax,nuvzmax,nx,ny,memind,
+ xt,yt,i,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 levinterpoln(uun,maxnests,nxmaxn,nymaxn,nuvzmax,ngrid,
+ nxn,nyn,memind,xtn,ytn,i,memtime(indexf),memtime(indexf+1),
+ itime,indexf,uint)
call levinterpoln(vvn,maxnests,nxmaxn,nymaxn,nuvzmax,ngrid,
+ nxn,nyn,memind,xtn,ytn,i,memtime(indexf),memtime(indexf+1),
+ itime,indexf,vint)
else ! linear interpolation
call levlininterpoln(uun,maxnests,nxmaxn,nymaxn,nuvzmax,ngrid,
+ nxn,nyn,memind,xtn,ytn,i,memtime(indexf),memtime(indexf+1),
+ itime,indexf,uint)
call levlininterpoln(vvn,maxnests,nxmaxn,nymaxn,nuvzmax,ngrid,
+ nxn,nyn,memind,xtn,ytn,i,memtime(indexf),memtime(indexf+1),
+ itime,indexf,vint)
endif
else ! polar stereographic grid
! only linear interpolation used
call levlininterpol(uupol,nxmax,nymax,nuvzmax,nx,ny,memind,
+ xt,yt,i,memtime(indexf),memtime(indexf+1),itime,indexf,
+ uint)
call levlininterpol(vvpol,nxmax,nymax,nuvzmax,nx,ny,memind,
+ xt,yt,i,memtime(indexf),memtime(indexf+1),itime,indexf,
+ vint)
endif
ph2=akm(i+1)+bkm(i+1)*psint
weight=weight+ph1-ph2
uhelp=uhelp+uint*(ph1-ph2)
vhelp=vhelp+vint*(ph1-ph2)
10 ph1=ph2
if (ngrid.ge.0) then ! mother grid
if ((inpolkind.eq.1).and.(iter.ne.1)) then
call levinterpol(uu,nxmax,nymax,nuvzmax,nx,ny,memind,
+ xt,yt,indz2,memtime(indexf),memtime(indexf+1),itime,indexf,
+ uint)
call levinterpol(vv,nxmax,nymax,nuvzmax,nx,ny,memind,
+ xt,yt,indz2,memtime(indexf),memtime(indexf+1),itime,indexf,
+ vint)
else ! bilinear interpolation
call levlininterpol(uu,nxmax,nymax,nuvzmax,nx,ny,memind,
+ xt,yt,indz2,memtime(indexf),memtime(indexf+1),itime,indexf,
+ uint)
call levlininterpol(vv,nxmax,nymax,nuvzmax,nx,ny,memind,
+ xt,yt,indz2,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 levinterpoln(uun,maxnests,nxmaxn,nymaxn,nuvzmax,ngrid,
+ nxn,nyn,memind,xtn,ytn,indz2,memtime(indexf),
+ memtime(indexf+1),itime,indexf,uint)
call levinterpoln(vvn,maxnests,nxmaxn,nymaxn,nuvzmax,ngrid,
+ nxn,nyn,memind,xtn,ytn,indz2,memtime(indexf),
+ memtime(indexf+1),itime,indexf,vint)
else ! linear interpolation
call levlininterpoln(uun,maxnests,nxmaxn,nymaxn,nuvzmax,ngrid,
+ nxn,nyn,memind,xtn,ytn,indz2,memtime(indexf),
+ memtime(indexf+1),itime,indexf,uint)
call levlininterpoln(vvn,maxnests,nxmaxn,nymaxn,nuvzmax,ngrid,
+ nxn,nyn,memind,xtn,ytn,indz2,memtime(indexf),
+ memtime(indexf+1),itime,indexf,vint)
endif
else ! polar stereographic grid
call levlininterpol(uupol,nxmax,nymax,nuvzmax,nx,ny,memind,
+ xt,yt,indz2,memtime(indexf),memtime(indexf+1),itime,indexf,
+ uint)
call levlininterpol(vvpol,nxmax,nymax,nuvzmax,nx,ny,memind,
+ xt,yt,indz2,memtime(indexf),memtime(indexf+1),itime,indexf,
+ vint)
endif
ph2=akz(indz2)+bkz(indz2)*psint
weight=weight+ph1-ph2
uhelp=uhelp+uint*(ph1-ph2)
vhelp=vhelp+vint*(ph1-ph2)
uint=uhelp/weight
vint=vhelp/weight
endif
return
end
|
