| 3 |
michaesp |
1 |
subroutine copyar(array1,array2,nsize)
|
|
|
2 |
C ======================================
|
|
|
3 |
C
|
|
|
4 |
C Copy the array1 on array2.
|
|
|
5 |
C
|
|
|
6 |
C array1 input array to copy
|
|
|
7 |
C array2 output array to write
|
|
|
8 |
C nsize input size of arrays
|
|
|
9 |
|
|
|
10 |
integer i,nsize
|
|
|
11 |
real array1(nsize),array2(nsize)
|
|
|
12 |
|
|
|
13 |
do i=1,nsize
|
|
|
14 |
array2(i)=array1(i)
|
|
|
15 |
enddo
|
|
|
16 |
|
|
|
17 |
return
|
|
|
18 |
end
|
|
|
19 |
subroutine clipreg2(plon,plat,varmin,varmax,nx,ny,ak,bk,nz,sp,
|
|
|
20 |
> prelevs,kcl0,kcl1,xcorner,ycorner,pmin,pmax,ind,
|
|
|
21 |
> nxy,single,onelev,ierr)
|
|
|
22 |
C ==============================================================
|
|
|
23 |
C
|
|
|
24 |
C Clip a region to calculate trajectories.
|
|
|
25 |
C
|
|
|
26 |
C plon input longitude of computational pole
|
|
|
27 |
C plat input latitude of computational pole
|
|
|
28 |
C varmin input array with minimum phys. coord. of data region
|
|
|
29 |
C varmax input array with maximum phys. coord. of data region
|
|
|
30 |
C nx input number of data points along latitude
|
|
|
31 |
C ny input number of data points along longitude
|
|
|
32 |
C ak input array contains ak values to calculate data-levels
|
|
|
33 |
C bk input array contains bk values to calculate data-levels
|
|
|
34 |
C nz input number of levels
|
|
|
35 |
C sp input array contains surface pressure
|
|
|
36 |
C prelevs input logical var (=.true. if data is on pressure levels)
|
|
|
37 |
C kcl0,1 output lowest,highest index for level of clipped region
|
|
|
38 |
C xcorner output x-coordinates of corners of clipped region
|
|
|
39 |
C ycorner output y-coordinates of corners of clipped region
|
|
|
40 |
C pmin output minimum pressure value of clipped region
|
|
|
41 |
C pmax output maximum pressure value of clipped region
|
|
|
42 |
C ind output array with index of grid points within clipped region
|
|
|
43 |
C nxy output number of points within clipped region
|
|
|
44 |
C single output logical flag for single trajectories
|
|
|
45 |
C onelev output logical flag if clipmin(3)=clipmax(3)
|
|
|
46 |
C ierr output error flag
|
|
|
47 |
|
|
|
48 |
real pmin,pmax,varmin(4),varmax(4)
|
|
|
49 |
real xcorner(4),ycorner(4),xmax,ymax,xmin,ymin,xx,yy
|
|
|
50 |
real ak(nz),bk(nz)
|
|
|
51 |
real sp(nx,ny)
|
|
|
52 |
real dx,dy
|
|
|
53 |
real plon,plat,xphys,yphys
|
|
|
54 |
integer kcl0,kcl1,i,j,k,l,nx,ny,nz,nxy,ierr
|
|
|
55 |
integer ind(400*100)
|
|
|
56 |
logical prelevs,onelev,single
|
|
|
57 |
|
|
|
58 |
real lmtolms,phtophs
|
|
|
59 |
|
|
|
60 |
C Reset error flag
|
|
|
61 |
|
|
|
62 |
ierr=0
|
|
|
63 |
|
|
|
64 |
C Calculate grid increments
|
|
|
65 |
|
|
|
66 |
dx=(varmax(1)-varmin(1))/(nx-1)
|
|
|
67 |
dy=(varmax(2)-varmin(2))/(ny-1)
|
|
|
68 |
|
|
|
69 |
C Read corners of region to clip from tape 9
|
|
|
70 |
|
|
|
71 |
read(9,10)xcorner(1)
|
|
|
72 |
read(9,10)ycorner(1)
|
|
|
73 |
read(9,10)xcorner(2)
|
|
|
74 |
read(9,10)ycorner(2)
|
|
|
75 |
read(9,10)xcorner(3)
|
|
|
76 |
read(9,10)ycorner(3)
|
|
|
77 |
read(9,10)xcorner(4)
|
|
|
78 |
read(9,10)ycorner(4)
|
|
|
79 |
read(9,10)pmax
|
|
|
80 |
read(9,10)pmin
|
|
|
81 |
|
|
|
82 |
10 format(f7.3)
|
|
|
83 |
|
|
|
84 |
C If necessary transform corners to the computational coordinates
|
|
|
85 |
|
|
|
86 |
if ((plon.ne.0.).or.(plat.ne.90.)) then
|
|
|
87 |
do i=1,4
|
|
|
88 |
xphys=lmtolms(ycorner(i),xcorner(i),plat,plon)
|
|
|
89 |
yphys=phtophs(ycorner(i),xcorner(i),plat,plon)
|
|
|
90 |
xcorner(i)=xphys
|
|
|
91 |
ycorner(i)=yphys
|
|
|
92 |
enddo
|
|
|
93 |
endif
|
|
|
94 |
|
|
|
95 |
C Test if corners are within data domain
|
|
|
96 |
|
|
|
97 |
do i=1,4
|
|
|
98 |
if (xcorner(i).lt.varmin(1)) goto 980
|
|
|
99 |
if (ycorner(i).lt.varmin(2)) goto 980
|
|
|
100 |
if (xcorner(i).gt.varmax(1)) goto 980
|
|
|
101 |
if (ycorner(i).gt.varmax(2)) goto 980
|
|
|
102 |
enddo
|
|
|
103 |
if (pmax.gt.varmin(3)) goto 980
|
|
|
104 |
if (pmin.lt.varmax(3)) goto 980
|
|
|
105 |
|
|
|
106 |
C Set onelev flag
|
|
|
107 |
|
|
|
108 |
if (pmin.eq.pmax) then
|
|
|
109 |
onelev=.true.
|
|
|
110 |
else
|
|
|
111 |
onelev=.false.
|
|
|
112 |
endif
|
|
|
113 |
|
|
|
114 |
C Check for single trajectory
|
|
|
115 |
|
|
|
116 |
single=.false.
|
|
|
117 |
if ((xcorner(1).eq.xcorner(2)).and.
|
|
|
118 |
> (xcorner(2).eq.xcorner(3)).and.
|
|
|
119 |
> (xcorner(3).eq.xcorner(4)).and.
|
|
|
120 |
> (ycorner(1).eq.ycorner(2)).and.
|
|
|
121 |
> (ycorner(2).eq.ycorner(3)).and.
|
|
|
122 |
> (ycorner(3).eq.ycorner(4))) then
|
|
|
123 |
if (pmin.eq.pmax) then
|
|
|
124 |
nxy=1
|
|
|
125 |
single=.true.
|
|
|
126 |
return
|
|
|
127 |
else
|
|
|
128 |
goto 980
|
|
|
129 |
endif
|
|
|
130 |
endif
|
|
|
131 |
|
|
|
132 |
C Determine which grid points are within clipped region
|
|
|
133 |
|
|
|
134 |
xmax=xcorner(1)
|
|
|
135 |
xmin=xcorner(1)
|
|
|
136 |
ymax=ycorner(1)
|
|
|
137 |
ymin=ycorner(1)
|
|
|
138 |
do i=2,4
|
|
|
139 |
if (xcorner(i).lt.xmin) xmin=xcorner(i)
|
|
|
140 |
if (xcorner(i).gt.xmax) xmax=xcorner(i)
|
|
|
141 |
if (ycorner(i).lt.ymin) ymin=ycorner(i)
|
|
|
142 |
if (ycorner(i).gt.ymax) ymax=ycorner(i)
|
|
|
143 |
enddo
|
|
|
144 |
|
|
|
145 |
nxy=0
|
|
|
146 |
|
|
|
147 |
do i=1,nx
|
|
|
148 |
xx=varmin(1)+(i-1)*dx
|
|
|
149 |
do j=1,ny
|
|
|
150 |
yy=varmin(2)+(j-1)*dy
|
|
|
151 |
if (xx.lt.xmin) goto 970
|
|
|
152 |
if (xx.gt.xmax) goto 970
|
|
|
153 |
if (yy.lt.ymin) goto 970
|
|
|
154 |
if (yy.gt.ymax) goto 970
|
|
|
155 |
if ((xx-xcorner(1))*(ycorner(2)-ycorner(1))-
|
|
|
156 |
> (yy-ycorner(1))*(xcorner(2)-xcorner(1)).gt.0.) goto 970
|
|
|
157 |
if ((xx-xcorner(2))*(ycorner(3)-ycorner(2))-
|
|
|
158 |
> (yy-ycorner(2))*(xcorner(3)-xcorner(2)).gt.0.) goto 970
|
|
|
159 |
if ((xx-xcorner(3))*(ycorner(4)-ycorner(3))-
|
|
|
160 |
> (yy-ycorner(3))*(xcorner(4)-xcorner(3)).gt.0.) goto 970
|
|
|
161 |
if ((xx-xcorner(4))*(ycorner(1)-ycorner(4))-
|
|
|
162 |
> (yy-ycorner(4))*(xcorner(1)-xcorner(4)).gt.0.) goto 970
|
|
|
163 |
|
|
|
164 |
nxy=nxy+1
|
|
|
165 |
ind(nxy)=i+(j-1)*nx
|
|
|
166 |
970 continue
|
|
|
167 |
enddo
|
|
|
168 |
enddo
|
|
|
169 |
|
|
|
170 |
C Calculate vertical indices for clipped region
|
|
|
171 |
|
|
|
172 |
if (.not.onelev) then
|
|
|
173 |
if (prelevs) then ! data on pressure levels
|
|
|
174 |
do k=1,nz-1
|
|
|
175 |
if ((pmax.le.ak(k)).and.(pmax.ge.ak(k+1))) then
|
|
|
176 |
if (ak(k)-pmax.gt.pmax-ak(k+1)) then
|
|
|
177 |
kcl0=k+1
|
|
|
178 |
else
|
|
|
179 |
kcl0=k
|
|
|
180 |
endif
|
|
|
181 |
endif
|
|
|
182 |
if ((pmin.le.ak(k)).and.(pmin.ge.ak(k+1))) then
|
|
|
183 |
if (ak(k)-pmin.gt.pmin-ak(k+1)) then
|
|
|
184 |
kcl1=k+1
|
|
|
185 |
else
|
|
|
186 |
kcl1=k
|
|
|
187 |
endif
|
|
|
188 |
endif
|
|
|
189 |
enddo
|
|
|
190 |
else ! data on model levels
|
|
|
191 |
kcl0=1 ! initialize kcl0 to lowest level
|
|
|
192 |
kcl1=nz ! initialize kcl1 to highest level
|
|
|
193 |
do k=nz-1,1,-1
|
|
|
194 |
do l=1,nxy
|
|
|
195 |
i=mod(mod(ind(l)-1,nx*ny),nx)+1
|
|
|
196 |
j=int(mod(ind(l)-1,nx*ny)/nx)+1
|
|
|
197 |
if ((pmax.le.ak(k)+bk(k)*sp(i,j)).and.
|
|
|
198 |
> (pmax.gt.ak(k+1)+bk(k+1)*sp(i,j))) kcl0=k+1
|
|
|
199 |
enddo
|
|
|
200 |
enddo
|
|
|
201 |
do k=1,nz-1
|
|
|
202 |
do l=1,nxy
|
|
|
203 |
i=mod(mod(ind(l)-1,nx*ny),nx)+1
|
|
|
204 |
j=int(mod(ind(l)-1,nx*ny)/nx)+1
|
|
|
205 |
if ((pmin.lt.ak(k)+bk(k)*sp(i,j)).and.
|
|
|
206 |
> (pmin.ge.ak(k+1)+bk(k+1)*sp(i,j))) kcl1=k
|
|
|
207 |
enddo
|
|
|
208 |
enddo
|
|
|
209 |
endif
|
|
|
210 |
endif
|
|
|
211 |
|
|
|
212 |
C Inform about eventually corrected boundary values
|
|
|
213 |
|
|
|
214 |
if (.not.onelev) then
|
|
|
215 |
if (prelevs) then
|
|
|
216 |
write(*,*)pmin,' corrected to --> ',ak(kcl0)
|
|
|
217 |
write(*,*)pmax,' corrected to --> ',ak(kcl1)
|
|
|
218 |
else
|
|
|
219 |
write(*,*)kcl0,' is level index for lower boundary'
|
|
|
220 |
write(*,*)kcl1,' is level index for upper boundary'
|
|
|
221 |
endif
|
|
|
222 |
endif
|
|
|
223 |
|
|
|
224 |
C Save eventually corrected boundary values
|
|
|
225 |
|
|
|
226 |
if (prelevs.and.(.not.onelev)) then
|
|
|
227 |
pmin=ak(kcl0)
|
|
|
228 |
pmax=ak(kcl1)
|
|
|
229 |
endif
|
|
|
230 |
|
|
|
231 |
return
|
|
|
232 |
|
|
|
233 |
980 ierr=1
|
|
|
234 |
return
|
|
|
235 |
end
|
|
|
236 |
subroutine getnpoints(plon,plat,varmin,varmax,nx,ny,ak,bk,nz,sp,
|
|
|
237 |
> levtyp,grid,delh,delp,vcflag,xcorner,ycorner,
|
|
|
238 |
> pmin,pmax,nxy,nlev,ierr)
|
|
|
239 |
C ================================================================
|
|
|
240 |
C
|
|
|
241 |
C Determine approximate number of points within clipped region.
|
|
|
242 |
C
|
|
|
243 |
C plon input longitude of computational pole
|
|
|
244 |
C plat input latitude of computational pole
|
|
|
245 |
C varmin input array with minimum phys. coord. of data region
|
|
|
246 |
C varmax input array with maximum phys. coord. of data region
|
|
|
247 |
C nx input number of data points along latitude
|
|
|
248 |
C ny input number of data points along longitude
|
|
|
249 |
C ak input array contains ak values to calculate data-levels
|
|
|
250 |
C bk input array contains bk values to calculate data-levels
|
|
|
251 |
C nz input number of levels
|
|
|
252 |
C sp input array contains surface pressure
|
|
|
253 |
C levtyp input int (=1 pressure =2 theta =3 model levels)
|
|
|
254 |
C delh input desired horizontal resol of starting points (deg or km)
|
|
|
255 |
C delp input desired vertical resol of starting points (hPa)
|
|
|
256 |
C xcorner output corners of clipped region
|
|
|
257 |
C pmin output minimum pressure value of clipped region
|
|
|
258 |
C pmax output maximum pressure value of clipped region
|
|
|
259 |
C nxy output number of points (horizontally) within clipped region
|
|
|
260 |
C nv output number of vertical points within clipped region
|
|
|
261 |
C ierr output error flag
|
|
|
262 |
|
|
|
263 |
real deltay,pi
|
|
|
264 |
parameter (deltay=111.2) ! distance in km between 2 lat circles
|
|
|
265 |
parameter (pi=3.1415927)
|
|
|
266 |
|
|
|
267 |
real pmin,pmax,varmin(4),varmax(4),delx,dely,delh,delp
|
|
|
268 |
real xcorner(4),ycorner(4),xmax,ymax,xmin,ymin,x,y
|
|
|
269 |
real ak(nz),bk(nz)
|
|
|
270 |
real sp(nx,ny)
|
|
|
271 |
real dx,dy
|
|
|
272 |
real plon,plat,xphys,yphys
|
|
|
273 |
integer kcl0,kcl1,i,j,k,nx,ny,nz,nxy,nlev,vcflag,ierr
|
|
|
274 |
integer nnx,nny
|
|
|
275 |
integer levtyp
|
|
|
276 |
logical grid
|
|
|
277 |
|
|
|
278 |
real lmtolms,phtophs
|
|
|
279 |
|
|
|
280 |
C Reset error flag
|
|
|
281 |
|
|
|
282 |
ierr=0
|
|
|
283 |
|
|
|
284 |
C Initialize nlev
|
|
|
285 |
|
|
|
286 |
nlev=0
|
|
|
287 |
|
|
|
288 |
C Calculate grid increments
|
|
|
289 |
|
|
|
290 |
dx=(varmax(1)-varmin(1))/(nx-1)
|
|
|
291 |
dy=(varmax(2)-varmin(2))/(ny-1)
|
|
|
292 |
|
|
|
293 |
C Read vcflag (=1 if pmin/max denote pressure values;
|
|
|
294 |
C =2 if pmin/max denote layers)
|
|
|
295 |
read(9,*)vcflag
|
|
|
296 |
|
|
|
297 |
C Read corners of region to clip from tape 9
|
|
|
298 |
|
|
|
299 |
read(9,10)xcorner(1)
|
|
|
300 |
read(9,10)ycorner(1)
|
|
|
301 |
read(9,10)xcorner(2)
|
|
|
302 |
read(9,10)ycorner(2)
|
|
|
303 |
read(9,10)xcorner(3)
|
|
|
304 |
read(9,10)ycorner(3)
|
|
|
305 |
read(9,10)xcorner(4)
|
|
|
306 |
read(9,10)ycorner(4)
|
|
|
307 |
read(9,10)pmax
|
|
|
308 |
read(9,10)pmin
|
|
|
309 |
10 format(f7.3)
|
|
|
310 |
|
|
|
311 |
if ((plon.eq.0.).and.(plat.eq.90.)) then
|
|
|
312 |
if (xcorner(1).eq.-999.) xcorner(1)=varmin(1)
|
|
|
313 |
if (xcorner(2).eq.-999.) xcorner(2)=varmax(1)-dx
|
|
|
314 |
if (xcorner(3).eq.-999.) xcorner(3)=varmax(1)-dx
|
|
|
315 |
if (xcorner(4).eq.-999.) xcorner(4)=varmin(1)
|
|
|
316 |
if (ycorner(1).eq.-999.) ycorner(1)=varmin(2)
|
|
|
317 |
if (ycorner(2).eq.-999.) ycorner(2)=varmin(2)
|
|
|
318 |
if (ycorner(3).eq.-999.) ycorner(3)=varmax(2)-dy
|
|
|
319 |
if (ycorner(4).eq.-999.) ycorner(4)=varmax(2)-dy
|
|
|
320 |
endif
|
|
|
321 |
|
|
|
322 |
C Special treatment for rotated coordinate systems (transform corners to
|
|
|
323 |
C rotated grid)
|
|
|
324 |
|
|
|
325 |
if ((plon.ne.0.).or.(abs(plat-90.).gt.1.e-3)) then
|
|
|
326 |
if (.not.grid) goto 975
|
|
|
327 |
do i=1,4
|
|
|
328 |
if ((xcorner(i).eq.-999.).and.(ycorner(i).eq.-999.)) then
|
|
|
329 |
if (i.eq.1) then
|
|
|
330 |
xcorner(i)=varmin(1)
|
|
|
331 |
ycorner(i)=varmin(2)
|
|
|
332 |
else if (i.eq.2) then
|
|
|
333 |
xcorner(i)=varmax(1)-dx
|
|
|
334 |
ycorner(i)=varmin(2)
|
|
|
335 |
else if (i.eq.3) then
|
|
|
336 |
xcorner(i)=varmax(1)-dx
|
|
|
337 |
ycorner(i)=varmax(2)-dy
|
|
|
338 |
else if (i.eq.4) then
|
|
|
339 |
xcorner(i)=varmin(1)
|
|
|
340 |
ycorner(i)=varmax(2)-dy
|
|
|
341 |
endif
|
|
|
342 |
else if ((xcorner(i).eq.-999.).or.(ycorner(i).eq.-999.)) then
|
|
|
343 |
goto 985
|
|
|
344 |
else
|
|
|
345 |
xphys=lmtolms(ycorner(i),xcorner(i),plat,plon)
|
|
|
346 |
yphys=phtophs(ycorner(i),xcorner(i),plat,plon)
|
|
|
347 |
xcorner(i)=xphys
|
|
|
348 |
ycorner(i)=yphys
|
|
|
349 |
endif
|
|
|
350 |
enddo
|
|
|
351 |
endif
|
|
|
352 |
|
|
|
353 |
C Test if corners are within data domain
|
|
|
354 |
|
|
|
355 |
do i=1,4
|
|
|
356 |
* if (xcorner(i).lt.varmin(1)) goto 980
|
|
|
357 |
* if (ycorner(i).lt.varmin(2)) goto 980
|
|
|
358 |
* if (xcorner(i).gt.varmax(1)) goto 980
|
|
|
359 |
* if (ycorner(i).gt.varmax(2)) goto 980
|
|
|
360 |
if (xcorner(i)+1.e-4.lt.varmin(1)) then
|
|
|
361 |
print*,'1 ',xcorner(i),varmin(1)
|
|
|
362 |
goto 980
|
|
|
363 |
endif
|
|
|
364 |
if (ycorner(i)+1.e-4.lt.varmin(2)) then
|
|
|
365 |
print*,'2 ',xcorner(i),varmin(2)
|
|
|
366 |
goto 980
|
|
|
367 |
endif
|
|
|
368 |
if (xcorner(i)-1.e-4.gt.varmax(1)) then
|
|
|
369 |
print*,'3 ',xcorner(i),varmax(1)
|
|
|
370 |
goto 980
|
|
|
371 |
endif
|
|
|
372 |
if (ycorner(i)-1.e-4.gt.varmax(2)) then
|
|
|
373 |
print*,'4 ',xcorner(i),varmax(2)
|
|
|
374 |
goto 980
|
|
|
375 |
endif
|
|
|
376 |
enddo
|
|
|
377 |
|
|
|
378 |
if (vcflag.eq.1) then
|
|
|
379 |
if (levtyp.eq.2) then
|
|
|
380 |
if (pmax.lt.varmin(3)) goto 980
|
|
|
381 |
if (pmin.gt.varmax(3)) goto 980
|
|
|
382 |
else
|
|
|
383 |
if (pmax.gt.varmin(3)) goto 980
|
|
|
384 |
if (pmin.lt.varmax(3)) goto 980
|
|
|
385 |
endif
|
|
|
386 |
endif
|
|
|
387 |
|
|
|
388 |
C Control output of clipped region
|
|
|
389 |
|
|
|
390 |
write(*,*)'corners of the clipped region:'
|
|
|
391 |
write(*,*)' ',xcorner(1),ycorner(1)
|
|
|
392 |
write(*,*)' ',xcorner(2),ycorner(2)
|
|
|
393 |
write(*,*)' ',xcorner(3),ycorner(3)
|
|
|
394 |
write(*,*)' ',xcorner(4),ycorner(4)
|
|
|
395 |
|
|
|
396 |
C Set nlev=1 if regions contains only one level
|
|
|
397 |
|
|
|
398 |
if (pmin.eq.pmax) nlev=1
|
|
|
399 |
|
|
|
400 |
C Check for single trajectory
|
|
|
401 |
|
|
|
402 |
if ((xcorner(1).eq.xcorner(2)).and.
|
|
|
403 |
> (xcorner(2).eq.xcorner(3)).and.
|
|
|
404 |
> (xcorner(3).eq.xcorner(4)).and.
|
|
|
405 |
> (ycorner(1).eq.ycorner(2)).and.
|
|
|
406 |
> (ycorner(2).eq.ycorner(3)).and.
|
|
|
407 |
> (ycorner(3).eq.ycorner(4))) then
|
|
|
408 |
if (nlev.eq.1) then
|
|
|
409 |
nxy=1
|
|
|
410 |
return
|
|
|
411 |
else
|
|
|
412 |
goto 980
|
|
|
413 |
endif
|
|
|
414 |
endif
|
|
|
415 |
|
|
|
416 |
C Determine number of grid points within clipped region (horizontally)
|
|
|
417 |
|
|
|
418 |
xmax=xcorner(1)
|
|
|
419 |
xmin=xcorner(1)
|
|
|
420 |
ymax=ycorner(1)
|
|
|
421 |
ymin=ycorner(1)
|
|
|
422 |
do i=2,4
|
|
|
423 |
if (xcorner(i).lt.xmin) xmin=xcorner(i)
|
|
|
424 |
if (xcorner(i).gt.xmax) xmax=xcorner(i)
|
|
|
425 |
if (ycorner(i).lt.ymin) ymin=ycorner(i)
|
|
|
426 |
if (ycorner(i).gt.ymax) ymax=ycorner(i)
|
|
|
427 |
enddo
|
|
|
428 |
|
|
|
429 |
nxy=0
|
|
|
430 |
|
|
|
431 |
if (grid) then
|
|
|
432 |
do i=1,nx
|
|
|
433 |
x=varmin(1)+(i-1)*dx
|
|
|
434 |
do j=1,ny
|
|
|
435 |
y=varmin(2)+(j-1)*dy
|
|
|
436 |
if (x.lt.xmin) goto 970
|
|
|
437 |
if (x.gt.xmax) goto 970
|
|
|
438 |
if (y.lt.ymin) goto 970
|
|
|
439 |
if (y.gt.ymax) goto 970
|
|
|
440 |
if ((x-xcorner(1))*(ycorner(2)-ycorner(1))-
|
|
|
441 |
> (y-ycorner(1))*(xcorner(2)-xcorner(1)).gt.0.) goto 970
|
|
|
442 |
if ((x-xcorner(2))*(ycorner(3)-ycorner(2))-
|
|
|
443 |
> (y-ycorner(2))*(xcorner(3)-xcorner(2)).gt.0.) goto 970
|
|
|
444 |
if ((x-xcorner(3))*(ycorner(4)-ycorner(3))-
|
|
|
445 |
> (y-ycorner(3))*(xcorner(4)-xcorner(3)).gt.0.) goto 970
|
|
|
446 |
if ((x-xcorner(4))*(ycorner(1)-ycorner(4))-
|
|
|
447 |
> (y-ycorner(4))*(xcorner(1)-xcorner(4)).gt.0.) goto 970
|
|
|
448 |
|
|
|
449 |
nxy=nxy+1
|
|
|
450 |
970 continue
|
|
|
451 |
enddo
|
|
|
452 |
enddo
|
|
|
453 |
else
|
|
|
454 |
nny=nint((ymax-ymin)*deltay/delh)+1
|
|
|
455 |
dely=(ymax-ymin)/real(nny-1)
|
|
|
456 |
print*,ymax,ymin,delh,nny
|
|
|
457 |
write(*,*)'value of dely set to ',dely,' (in deg)'
|
|
|
458 |
do j=1,nny
|
|
|
459 |
y=ymin+(j-1)*dely
|
|
|
460 |
nnx=nint((xmax-xmin)*cos(y*pi/180.)*deltay/delh)
|
|
|
461 |
if (nnx.gt.0) then
|
|
|
462 |
delx=(xmax-xmin)/real(nnx)
|
|
|
463 |
else
|
|
|
464 |
delx=360.
|
|
|
465 |
nnx=1
|
|
|
466 |
endif
|
|
|
467 |
write(*,*)'val of delx at lat ',y,' set to ',delx,' (in deg)'
|
|
|
468 |
do i=1,nnx
|
|
|
469 |
x=xmin+(i-1)*delx
|
|
|
470 |
if ((x-xcorner(1))*(ycorner(2)-ycorner(1))-
|
|
|
471 |
> (y-ycorner(1))*(xcorner(2)-xcorner(1)).gt.0.01) goto 971
|
|
|
472 |
if ((x-xcorner(2))*(ycorner(3)-ycorner(2))-
|
|
|
473 |
> (y-ycorner(2))*(xcorner(3)-xcorner(2)).gt.0.01) goto 971
|
|
|
474 |
if ((x-xcorner(3))*(ycorner(4)-ycorner(3))-
|
|
|
475 |
> (y-ycorner(3))*(xcorner(4)-xcorner(3)).gt.0.01) goto 971
|
|
|
476 |
if ((x-xcorner(4))*(ycorner(1)-ycorner(4))-
|
|
|
477 |
> (y-ycorner(4))*(xcorner(1)-xcorner(4)).gt.0.01) goto 971
|
|
|
478 |
|
|
|
479 |
nxy=nxy+1
|
|
|
480 |
971 continue
|
|
|
481 |
enddo
|
|
|
482 |
enddo
|
|
|
483 |
endif
|
|
|
484 |
|
|
|
485 |
write(*,*)'num of pts in the clipped region (horizontally): ',nxy
|
|
|
486 |
|
|
|
487 |
if (.not.grid) then
|
|
|
488 |
nlev=nint((pmax-pmin)/delp)+1
|
|
|
489 |
write(*,*)'number of levels in the data domain: ',nlev
|
|
|
490 |
return
|
|
|
491 |
endif
|
|
|
492 |
|
|
|
493 |
C Calculate vertical indices for whole data domain
|
|
|
494 |
|
|
|
495 |
if (vcflag.eq.1) then
|
|
|
496 |
if (nlev.ne.1) then
|
|
|
497 |
if (levtyp.eq.1) then ! data on pressure levels
|
|
|
498 |
do k=1,nz-1
|
|
|
499 |
if ((pmax.le.ak(k)).and.(pmax.ge.ak(k+1))) then
|
|
|
500 |
if (ak(k)-pmax.gt.pmax-ak(k+1)) then
|
|
|
501 |
kcl0=k+1
|
|
|
502 |
else
|
|
|
503 |
kcl0=k
|
|
|
504 |
endif
|
|
|
505 |
endif
|
|
|
506 |
if ((pmin.le.ak(k)).and.(pmin.ge.ak(k+1))) then
|
|
|
507 |
if (ak(k)-pmin.gt.pmin-ak(k+1)) then
|
|
|
508 |
kcl1=k+1
|
|
|
509 |
else
|
|
|
510 |
kcl1=k
|
|
|
511 |
endif
|
|
|
512 |
endif
|
|
|
513 |
enddo
|
|
|
514 |
else if (levtyp.eq.2) then ! data on theta levels
|
|
|
515 |
do k=1,nz-1
|
|
|
516 |
if ((pmax.ge.ak(k)).and.(pmax.le.ak(k+1))) then
|
|
|
517 |
if (pmax-ak(k).gt.ak(k+1)-pmax) then
|
|
|
518 |
kcl0=k+1
|
|
|
519 |
else
|
|
|
520 |
kcl0=k
|
|
|
521 |
endif
|
|
|
522 |
endif
|
|
|
523 |
if ((pmin.ge.ak(k)).and.(pmin.le.ak(k+1))) then
|
|
|
524 |
if (pmin-ak(k).gt.ak(k+1)-pmin) then
|
|
|
525 |
kcl1=k+1
|
|
|
526 |
else
|
|
|
527 |
kcl1=k
|
|
|
528 |
endif
|
|
|
529 |
endif
|
|
|
530 |
enddo
|
|
|
531 |
else ! data on model levels
|
|
|
532 |
kcl0=1 ! initialize kcl0 to lowest level
|
|
|
533 |
kcl1=nz ! initialize kcl1 to highest level
|
|
|
534 |
do k=nz-1,1,-1
|
|
|
535 |
do i=1,nx
|
|
|
536 |
do j=1,ny
|
|
|
537 |
if ((pmax.le.ak(k)+bk(k)*sp(i,j)).and.
|
|
|
538 |
> (pmax.gt.ak(k+1)+bk(k+1)*sp(i,j))) then
|
|
|
539 |
kcl0=k+1
|
|
|
540 |
goto 950
|
|
|
541 |
endif
|
|
|
542 |
enddo
|
|
|
543 |
enddo
|
|
|
544 |
950 continue
|
|
|
545 |
enddo
|
|
|
546 |
if (kcl0.eq.2) then
|
|
|
547 |
do i=1,nx
|
|
|
548 |
do j=1,ny
|
|
|
549 |
if (pmax.ge.ak(1)+bk(1)*sp(i,j)) then
|
|
|
550 |
kcl0=1
|
|
|
551 |
goto 952
|
|
|
552 |
endif
|
|
|
553 |
enddo
|
|
|
554 |
enddo
|
|
|
555 |
endif
|
|
|
556 |
952 continue
|
|
|
557 |
do k=1,nz-1
|
|
|
558 |
do i=1,nx
|
|
|
559 |
do j=1,ny
|
|
|
560 |
if ((pmin.lt.ak(k)+bk(k)*sp(i,j)).and.
|
|
|
561 |
> (pmin.ge.ak(k+1)+bk(k+1)*sp(i,j))) then
|
|
|
562 |
kcl1=k
|
|
|
563 |
goto 951
|
|
|
564 |
endif
|
|
|
565 |
enddo
|
|
|
566 |
enddo
|
|
|
567 |
951 continue
|
|
|
568 |
enddo
|
|
|
569 |
endif
|
|
|
570 |
nlev=kcl1-kcl0+1
|
|
|
571 |
write(*,*)'number of levels in the data domain: ',nlev
|
|
|
572 |
write(*,*)'from indices',kcl0,' to ',kcl1
|
|
|
573 |
endif
|
|
|
574 |
else if (vcflag.eq.2) then
|
|
|
575 |
kcl0=nint(pmin)
|
|
|
576 |
kcl1=nint(pmax)
|
|
|
577 |
nlev=kcl1-kcl0+1
|
|
|
578 |
write(*,*)'number of levels in the data domain: ',nlev
|
|
|
579 |
write(*,*)'from indices',kcl0,' to ',kcl1
|
|
|
580 |
endif
|
|
|
581 |
|
|
|
582 |
return
|
|
|
583 |
|
|
|
584 |
975 ierr=3
|
|
|
585 |
return
|
|
|
586 |
980 ierr=1
|
|
|
587 |
return
|
|
|
588 |
985 ierr=2
|
|
|
589 |
return
|
|
|
590 |
end
|
|
|
591 |
subroutine clipreg(plon,plat,varmin,varmax,nx,ny,ak,bk,nz,sp,
|
|
|
592 |
> levtyp,grid,delh,delp,vcflag,xcorner,ycorner,
|
|
|
593 |
> pmin,pmax,xx,yy,pp,ntra,ierr)
|
|
|
594 |
C =============================================================
|
|
|
595 |
C
|
|
|
596 |
C Clip a region to calculate trajectories.
|
|
|
597 |
C
|
|
|
598 |
C plon input longitude of computational pole
|
|
|
599 |
C plat input latitude of computational pole
|
|
|
600 |
C varmin input array with minimum phys. coord. of data region
|
|
|
601 |
C varmax input array with maximum phys. coord. of data region
|
|
|
602 |
C nx input number of data points along latitude
|
|
|
603 |
C ny input number of data points along longitude
|
|
|
604 |
C ak input array contains ak values to calculate data-levels
|
|
|
605 |
C bk input array contains bk values to calculate data-levels
|
|
|
606 |
C nz input number of levels
|
|
|
607 |
C sp input array contains surface pressure
|
|
|
608 |
C levtyp input int (=1 pressure =2 theta =3 model levels)
|
|
|
609 |
C delh input desired horizontal resol of starting points (km)
|
|
|
610 |
C delp input desired vertical resol of starting points (hPa)
|
|
|
611 |
C xx output starting coordinate of trajectories (longitude)
|
|
|
612 |
C yy output starting coordinate of trajectories (latitude)
|
|
|
613 |
C pp output starting coordinate of trajectories (pressure)
|
|
|
614 |
C ntra in/output number of grid points within clipped region
|
|
|
615 |
C ierr output error flag
|
|
|
616 |
|
|
|
617 |
real deltay,pi
|
|
|
618 |
parameter (deltay=111.2) ! distance in km between 2 lat circles
|
|
|
619 |
parameter (pi=3.1415927)
|
|
|
620 |
|
|
|
621 |
real,allocatable, dimension (:) :: xval,yval,spval
|
|
|
622 |
|
|
|
623 |
real xx(*),yy(*),pp(*)
|
|
|
624 |
real pmin,pmax,pval,varmin(4),varmax(4),delh,delx,dely,delp
|
|
|
625 |
real xcorner(4),ycorner(4),xmax,ymax,xmin,ymin,x,y
|
|
|
626 |
real ak(nz),bk(nz)
|
|
|
627 |
real sp(nx,ny)
|
|
|
628 |
real dx,dy
|
|
|
629 |
real plon,plat,xphys,yphys
|
|
|
630 |
integer kcl0,kcl1,i,j,k,l,nx,ny,nz,nxy,nlev,ntra,vcflag,ierr
|
|
|
631 |
integer ii,jj
|
|
|
632 |
integer levtyp,stat
|
|
|
633 |
logical grid
|
|
|
634 |
|
|
|
635 |
real lmtolms,phtophs,int2d
|
|
|
636 |
|
|
|
637 |
C Allocate memory for dynamical arrays
|
|
|
638 |
|
|
|
639 |
allocate(xval(ntra),stat=stat)
|
|
|
640 |
if (stat.ne.0) print*,'*** error allocating array xval ***'
|
|
|
641 |
allocate(yval(ntra),stat=stat)
|
|
|
642 |
if (stat.ne.0) print*,'*** error allocating array yval ***'
|
|
|
643 |
allocate(spval(ntra),stat=stat)
|
|
|
644 |
if (stat.ne.0) print*,'*** error allocating array spval ***'
|
|
|
645 |
|
|
|
646 |
C Reset error flag
|
|
|
647 |
|
|
|
648 |
ierr=0
|
|
|
649 |
|
|
|
650 |
C Calculate grid increments
|
|
|
651 |
|
|
|
652 |
dx=(varmax(1)-varmin(1))/(nx-1)
|
|
|
653 |
dy=(varmax(2)-varmin(2))/(ny-1)
|
|
|
654 |
|
|
|
655 |
write(*,*)'corners of the clipped region in SR clipreg:'
|
|
|
656 |
write(*,*)' ',xcorner(1),ycorner(1)
|
|
|
657 |
write(*,*)' ',xcorner(2),ycorner(2)
|
|
|
658 |
write(*,*)' ',xcorner(3),ycorner(3)
|
|
|
659 |
write(*,*)' ',xcorner(4),ycorner(4)
|
|
|
660 |
|
|
|
661 |
C Set nlev=1 if regions contains only one level
|
|
|
662 |
|
|
|
663 |
if (pmin.eq.pmax) nlev=1
|
|
|
664 |
|
|
|
665 |
C Check for single trajectory
|
|
|
666 |
|
|
|
667 |
if ((xcorner(1).eq.xcorner(2)).and.
|
|
|
668 |
> (xcorner(2).eq.xcorner(3)).and.
|
|
|
669 |
> (xcorner(3).eq.xcorner(4)).and.
|
|
|
670 |
> (ycorner(1).eq.ycorner(2)).and.
|
|
|
671 |
> (ycorner(2).eq.ycorner(3)).and.
|
|
|
672 |
> (ycorner(3).eq.ycorner(4))) then
|
|
|
673 |
if (nlev.eq.1) then
|
|
|
674 |
ntra=1
|
|
|
675 |
xx(1)=xcorner(1)
|
|
|
676 |
yy(1)=ycorner(1)
|
|
|
677 |
if (vcflag.eq.1) then
|
|
|
678 |
pp(1)=pmin
|
|
|
679 |
else if (vcflag.eq.2) then
|
|
|
680 |
ii=nint((xx(1)-varmin(1))/dx)+1
|
|
|
681 |
jj=nint((yy(1)-varmin(2))/dy)+1
|
|
|
682 |
pp(1)=ak(nint(pmin))+bk(nint(pmin))*sp(ii,jj)
|
|
|
683 |
endif
|
|
|
684 |
return
|
|
|
685 |
else
|
|
|
686 |
goto 980
|
|
|
687 |
endif
|
|
|
688 |
endif
|
|
|
689 |
|
|
|
690 |
C Determine which grid points are within clipped region (horizontally)
|
|
|
691 |
|
|
|
692 |
xmax=xcorner(1)
|
|
|
693 |
xmin=xcorner(1)
|
|
|
694 |
ymax=ycorner(1)
|
|
|
695 |
ymin=ycorner(1)
|
|
|
696 |
do i=2,4
|
|
|
697 |
if (xcorner(i).lt.xmin) xmin=xcorner(i)
|
|
|
698 |
if (xcorner(i).gt.xmax) xmax=xcorner(i)
|
|
|
699 |
if (ycorner(i).lt.ymin) ymin=ycorner(i)
|
|
|
700 |
if (ycorner(i).gt.ymax) ymax=ycorner(i)
|
|
|
701 |
enddo
|
|
|
702 |
|
|
|
703 |
C Cut the last grid points of the data domain
|
|
|
704 |
C (this prevents problems in subroutine trace: calls to getsdat)
|
|
|
705 |
|
|
|
706 |
* if (xmax.gt.varmax(1)-dx) xmax=varmax(1)-dx
|
|
|
707 |
* if (ymax.gt.varmax(2)-dy) ymax=varmax(2)-dy
|
|
|
708 |
|
|
|
709 |
nxy=0
|
|
|
710 |
|
|
|
711 |
if (grid) then
|
|
|
712 |
do i=1,nx
|
|
|
713 |
x=varmin(1)+(i-1)*dx
|
|
|
714 |
do j=1,ny
|
|
|
715 |
y=varmin(2)+(j-1)*dy
|
|
|
716 |
* print*,x,y,xmin,xmax,ymin,ymax,
|
|
|
717 |
* > xcorner(1),ycorner(1),xcorner(2),ycorner(2),
|
|
|
718 |
* > xcorner(3),ycorner(3),xcorner(4),ycorner(4)
|
|
|
719 |
if (x.lt.xmin) goto 970
|
|
|
720 |
if (x.gt.xmax) goto 970
|
|
|
721 |
if (y.lt.ymin) goto 970
|
|
|
722 |
if (y.gt.ymax) goto 970
|
|
|
723 |
if ((x-xcorner(1))*(ycorner(2)-ycorner(1))-
|
|
|
724 |
> (y-ycorner(1))*(xcorner(2)-xcorner(1)).gt.0.) goto 970
|
|
|
725 |
if ((x-xcorner(2))*(ycorner(3)-ycorner(2))-
|
|
|
726 |
> (y-ycorner(2))*(xcorner(3)-xcorner(2)).gt.0.) goto 970
|
|
|
727 |
if ((x-xcorner(3))*(ycorner(4)-ycorner(3))-
|
|
|
728 |
> (y-ycorner(3))*(xcorner(4)-xcorner(3)).gt.0.) goto 970
|
|
|
729 |
if ((x-xcorner(4))*(ycorner(1)-ycorner(4))-
|
|
|
730 |
> (y-ycorner(4))*(xcorner(1)-xcorner(4)).gt.0.) goto 970
|
|
|
731 |
|
|
|
732 |
nxy=nxy+1
|
|
|
733 |
* print*,i,j,nxy
|
|
|
734 |
xval(nxy)=x
|
|
|
735 |
yval(nxy)=y
|
|
|
736 |
spval(nxy)=sp(i,j)
|
|
|
737 |
970 continue
|
|
|
738 |
enddo
|
|
|
739 |
enddo
|
|
|
740 |
else
|
|
|
741 |
nny=nint((ymax-ymin)*deltay/delh)+1
|
|
|
742 |
dely=(ymax-ymin)/real(nny-1)
|
|
|
743 |
print*,ymax,ymin,delh,nny
|
|
|
744 |
write(*,*)'value of dely set to ',dely,' (in deg)'
|
|
|
745 |
do j=1,nny
|
|
|
746 |
y=ymin+(j-1)*dely
|
|
|
747 |
nnx=nint((xmax-xmin)*cos(y*pi/180.)*deltay/delh)
|
|
|
748 |
if (nnx.gt.0) then
|
|
|
749 |
delx=(xmax-xmin)/real(nnx)
|
|
|
750 |
else
|
|
|
751 |
delx=360.
|
|
|
752 |
nnx=1
|
|
|
753 |
endif
|
|
|
754 |
write(*,*)'val of delx at lat ',y,' set to ',delx,' (in deg)'
|
|
|
755 |
do i=1,nnx
|
|
|
756 |
x=xmin+(i-1)*delx
|
|
|
757 |
if ((x-xcorner(1))*(ycorner(2)-ycorner(1))-
|
|
|
758 |
> (y-ycorner(1))*(xcorner(2)-xcorner(1)).gt.0.01) goto 971
|
|
|
759 |
if ((x-xcorner(2))*(ycorner(3)-ycorner(2))-
|
|
|
760 |
> (y-ycorner(2))*(xcorner(3)-xcorner(2)).gt.0.01) goto 971
|
|
|
761 |
if ((x-xcorner(3))*(ycorner(4)-ycorner(3))-
|
|
|
762 |
> (y-ycorner(3))*(xcorner(4)-xcorner(3)).gt.0.01) goto 971
|
|
|
763 |
if ((x-xcorner(4))*(ycorner(1)-ycorner(4))-
|
|
|
764 |
> (y-ycorner(4))*(xcorner(1)-xcorner(4)).gt.0.01) goto 971
|
|
|
765 |
|
|
|
766 |
nxy=nxy+1
|
|
|
767 |
* print*,i,j,nxy
|
|
|
768 |
xval(nxy)=x
|
|
|
769 |
yval(nxy)=y
|
|
|
770 |
* write(*,'(i8,2f9.2)')nxy,xval(nxy),yval(nxy)
|
|
|
771 |
ri=(x-varmin(1))/dx+1.
|
|
|
772 |
rj=(y-varmin(2))/dy+1.
|
|
|
773 |
spval(nxy)=int2d(sp,nx,ny,ri,rj)
|
|
|
774 |
971 continue
|
|
|
775 |
enddo
|
|
|
776 |
enddo
|
|
|
777 |
endif
|
|
|
778 |
write(*,*)'num of pts in the clipped region (horizontally): ',nxy
|
|
|
779 |
|
|
|
780 |
if (.not.grid) then
|
|
|
781 |
nlev=nint((pmax-pmin)/delp)+1
|
|
|
782 |
delp=(pmax-pmin)/real(nlev-1)
|
|
|
783 |
write(*,*)'number of levels in the data domain: ',nlev
|
|
|
784 |
goto 960
|
|
|
785 |
endif
|
|
|
786 |
|
|
|
787 |
C Calculate vertical indices for clipped region
|
|
|
788 |
|
|
|
789 |
if (vcflag.eq.1) then
|
|
|
790 |
if (nlev.ne.1) then
|
|
|
791 |
if (levtyp.eq.1) then ! data on pressure levels
|
|
|
792 |
do k=1,nz-1
|
|
|
793 |
if ((pmax.le.ak(k)).and.(pmax.ge.ak(k+1))) then
|
|
|
794 |
if (ak(k)-pmax.gt.pmax-ak(k+1)) then
|
|
|
795 |
kcl0=k+1
|
|
|
796 |
else
|
|
|
797 |
kcl0=k
|
|
|
798 |
endif
|
|
|
799 |
endif
|
|
|
800 |
if ((pmin.le.ak(k)).and.(pmin.ge.ak(k+1))) then
|
|
|
801 |
if (ak(k)-pmin.gt.pmin-ak(k+1)) then
|
|
|
802 |
kcl1=k+1
|
|
|
803 |
else
|
|
|
804 |
kcl1=k
|
|
|
805 |
endif
|
|
|
806 |
endif
|
|
|
807 |
enddo
|
|
|
808 |
else if (levtyp.eq.2) then ! data on theta levels
|
|
|
809 |
do k=1,nz-1
|
|
|
810 |
if ((pmax.ge.ak(k)).and.(pmax.le.ak(k+1))) then
|
|
|
811 |
if (pmax-ak(k).gt.ak(k+1)-pmax) then
|
|
|
812 |
kcl0=k+1
|
|
|
813 |
else
|
|
|
814 |
kcl0=k
|
|
|
815 |
endif
|
|
|
816 |
endif
|
|
|
817 |
if ((pmin.ge.ak(k)).and.(pmin.le.ak(k+1))) then
|
|
|
818 |
if (pmin-ak(k).gt.ak(k+1)-pmin) then
|
|
|
819 |
kcl1=k+1
|
|
|
820 |
else
|
|
|
821 |
kcl1=k
|
|
|
822 |
endif
|
|
|
823 |
endif
|
|
|
824 |
enddo
|
|
|
825 |
else ! data on model levels
|
|
|
826 |
kcl0=1 ! initialize kcl0 to lowest level
|
|
|
827 |
kcl1=nz ! initialize kcl1 to highest level
|
|
|
828 |
do k=nz-1,1,-1
|
|
|
829 |
do l=1,nxy
|
|
|
830 |
if ((pmax.le.ak(k)+bk(k)*spval(l)).and.
|
|
|
831 |
> (pmax.gt.ak(k+1)+bk(k+1)*spval(l))) then
|
|
|
832 |
kcl0=k+1
|
|
|
833 |
goto 950
|
|
|
834 |
endif
|
|
|
835 |
enddo
|
|
|
836 |
950 continue
|
|
|
837 |
enddo
|
|
|
838 |
if (kcl0.eq.2) then
|
|
|
839 |
do l=1,nxy
|
|
|
840 |
if (pmax.ge.ak(1)+bk(1)*spval(l)) then
|
|
|
841 |
kcl0=1
|
|
|
842 |
goto 952
|
|
|
843 |
endif
|
|
|
844 |
enddo
|
|
|
845 |
endif
|
|
|
846 |
952 continue
|
|
|
847 |
do k=1,nz-1
|
|
|
848 |
do l=1,nxy
|
|
|
849 |
if ((pmin.lt.ak(k)+bk(k)*spval(l)).and.
|
|
|
850 |
> (pmin.ge.ak(k+1)+bk(k+1)*spval(l))) then
|
|
|
851 |
kcl1=k
|
|
|
852 |
goto 951
|
|
|
853 |
endif
|
|
|
854 |
enddo
|
|
|
855 |
951 continue
|
|
|
856 |
enddo
|
|
|
857 |
endif
|
|
|
858 |
nlev=kcl1-kcl0+1
|
|
|
859 |
write(*,*)'number of levels in the clipped region: ',nlev
|
|
|
860 |
write(*,*)'from indices ',kcl0,' to ',kcl1
|
|
|
861 |
endif
|
|
|
862 |
else if (vcflag.eq.2) then
|
|
|
863 |
kcl0=nint(pmin)
|
|
|
864 |
kcl1=nint(pmax)
|
|
|
865 |
nlev=kcl1-kcl0+1
|
|
|
866 |
write(*,*)'number of levels in the clipped region: ',nlev
|
|
|
867 |
write(*,*)'from indices ',kcl0,' to ',kcl1
|
|
|
868 |
endif
|
|
|
869 |
|
|
|
870 |
960 continue
|
|
|
871 |
|
|
|
872 |
C Serach grid points within clipped region
|
|
|
873 |
|
|
|
874 |
ntra=0
|
|
|
875 |
|
|
|
876 |
do k=1,nlev
|
|
|
877 |
do l=1,nxy
|
|
|
878 |
if (nlev.eq.1) then
|
|
|
879 |
pval=pmin
|
|
|
880 |
else
|
|
|
881 |
if (grid) then
|
|
|
882 |
pval=ak(k+kcl0-1)+bk(k+kcl0-1)*spval(l)
|
|
|
883 |
else
|
|
|
884 |
pval=pmin+(k-1)*delp
|
|
|
885 |
endif
|
|
|
886 |
endif
|
|
|
887 |
|
|
|
888 |
if (vcflag.eq.1) then
|
|
|
889 |
if ((levtyp.eq.1).and.(pval.le.spval(l))) then
|
|
|
890 |
ntra=ntra+1
|
|
|
891 |
xx(ntra)=xval(l)
|
|
|
892 |
yy(ntra)=yval(l)
|
|
|
893 |
pp(ntra)=pval
|
|
|
894 |
else if (levtyp.eq.2) then
|
|
|
895 |
ntra=ntra+1
|
|
|
896 |
xx(ntra)=xval(l)
|
|
|
897 |
yy(ntra)=yval(l)
|
|
|
898 |
pp(ntra)=pval
|
|
|
899 |
else if (levtyp.eq.3) then
|
|
|
900 |
if (grid) then
|
|
|
901 |
if ((pval.ge.pmin).and.(pval.le.pmax)) then
|
|
|
902 |
ntra=ntra+1
|
|
|
903 |
xx(ntra)=xval(l)
|
|
|
904 |
yy(ntra)=yval(l)
|
|
|
905 |
pp(ntra)=pval
|
|
|
906 |
endif
|
|
|
907 |
else
|
|
|
908 |
if (pval.le.spval(l)) then
|
|
|
909 |
ntra=ntra+1
|
|
|
910 |
xx(ntra)=xval(l)
|
|
|
911 |
yy(ntra)=yval(l)
|
|
|
912 |
pp(ntra)=pval
|
|
|
913 |
endif
|
|
|
914 |
endif
|
|
|
915 |
endif
|
|
|
916 |
else if (vcflag.eq.2) then
|
|
|
917 |
if ((levtyp.eq.1).and.(pval.le.spval(l))) then
|
|
|
918 |
ntra=ntra+1
|
|
|
919 |
xx(ntra)=xval(l)
|
|
|
920 |
yy(ntra)=yval(l)
|
|
|
921 |
pp(ntra)=pval
|
|
|
922 |
else
|
|
|
923 |
ntra=ntra+1
|
|
|
924 |
xx(ntra)=xval(l)
|
|
|
925 |
yy(ntra)=yval(l)
|
|
|
926 |
pp(ntra)=pval
|
|
|
927 |
endif
|
|
|
928 |
endif
|
|
|
929 |
enddo
|
|
|
930 |
enddo
|
|
|
931 |
|
|
|
932 |
write(*,*)'number of grid points in clipped region: ',ntra
|
|
|
933 |
|
|
|
934 |
C Inform about eventually corrected boundary values
|
|
|
935 |
|
|
|
936 |
if ((levtyp.eq.1).and.(nlev.gt.1)) then
|
|
|
937 |
write(*,*)pmin,' corrected from ',pmin,' to --> ',ak(kcl0)
|
|
|
938 |
write(*,*)pmax,' corrected from ',pmax,' to --> ',ak(kcl1)
|
|
|
939 |
pmin=ak(kcl0)
|
|
|
940 |
pmax=ak(kcl1)
|
|
|
941 |
endif
|
|
|
942 |
|
|
|
943 |
return
|
|
|
944 |
|
|
|
945 |
980 ierr=1
|
|
|
946 |
return
|
|
|
947 |
end
|
|
|
948 |
subroutine timediff(date1,date2,diff)
|
|
|
949 |
C =====================================
|
|
|
950 |
|
|
|
951 |
C New version with hour and minutes! (for hour and step [in hours]
|
|
|
952 |
C use the routine oldtimediff!)
|
|
|
953 |
C
|
|
|
954 |
C Calculates the time difference in hours (and minutes) for the two
|
|
|
955 |
C dates specified by the two arrays date1 and date2.
|
|
|
956 |
C They are expected to contain the following date information:
|
|
|
957 |
C year month day hour minute.
|
|
|
958 |
C
|
|
|
959 |
C date1 array specifying the first date
|
|
|
960 |
C date2 array specifying the second date
|
|
|
961 |
C diff time differenc between date1 and date2 in hours
|
|
|
962 |
C
|
|
|
963 |
|
|
|
964 |
integer date1(5),date2(5)
|
|
|
965 |
integer idays(12) ! array containing the days of the monthes
|
|
|
966 |
real diff
|
|
|
967 |
integer ixday,imdiff,ihdiff,iddiff,j
|
|
|
968 |
integer yy,yy1,yy2
|
|
|
969 |
|
|
|
970 |
idays(1)=31
|
|
|
971 |
idays(2)=28
|
|
|
972 |
idays(3)=31
|
|
|
973 |
idays(4)=30
|
|
|
974 |
idays(5)=31
|
|
|
975 |
idays(6)=30
|
|
|
976 |
idays(7)=31
|
|
|
977 |
idays(8)=31
|
|
|
978 |
idays(9)=30
|
|
|
979 |
idays(10)=31
|
|
|
980 |
idays(11)=30
|
|
|
981 |
idays(12)=31
|
|
|
982 |
|
|
|
983 |
C Check format of year (YYYY or YY - in case of YY assume 19YY)
|
|
|
984 |
|
|
|
985 |
if (date1(1).lt.100) date1(1)=1900+date1(1)
|
|
|
986 |
if (date2(1).lt.100) date2(1)=1900+date2(1)
|
|
|
987 |
|
|
|
988 |
C Determine if the period between date1 and date2 contains a Feb.29
|
|
|
989 |
|
|
|
990 |
ixday=0 ! extra day flag
|
|
|
991 |
|
|
|
992 |
yy1=min(date1(1),date2(1))
|
|
|
993 |
yy2=max(date1(1),date2(1))
|
|
|
994 |
if (yy1.eq.yy2) then
|
|
|
995 |
if (mod(yy1,4).eq.0) then
|
|
|
996 |
idays(2)=29
|
|
|
997 |
endif
|
|
|
998 |
else
|
|
|
999 |
if (mod(yy1,4).eq.0) then
|
|
|
1000 |
if (((yy1.eq.date1(1)).and.(date1(2).le.2)).or.
|
|
|
1001 |
> ((yy1.eq.date2(1)).and.(date2(2).le.2))) then
|
|
|
1002 |
ixday=ixday+1
|
|
|
1003 |
endif
|
|
|
1004 |
endif
|
|
|
1005 |
if (mod(yy2,4).eq.0) then
|
|
|
1006 |
if (((yy2.eq.date1(1)).and.(date1(2).gt.2)).or.
|
|
|
1007 |
> ((yy2.eq.date2(1)).and.(date2(2).gt.2))) then
|
|
|
1008 |
ixday=ixday+1
|
|
|
1009 |
endif
|
|
|
1010 |
endif
|
|
|
1011 |
if (yy2-yy1.gt.1) then
|
|
|
1012 |
do yy=yy1+1,yy2-1
|
|
|
1013 |
if (mod(yy,4).eq.0) then
|
|
|
1014 |
ixday=ixday+1
|
|
|
1015 |
endif
|
|
|
1016 |
enddo
|
|
|
1017 |
endif
|
|
|
1018 |
endif
|
|
|
1019 |
|
|
|
1020 |
ihdiff=0 ! diff. in hours between date1/date2
|
|
|
1021 |
iddiff=0 ! diff. in days between date1/date2
|
|
|
1022 |
|
|
|
1023 |
if (date1(1).gt.date2(1)) then ! compare years
|
|
|
1024 |
do j=date2(1),date1(1)-1
|
|
|
1025 |
iddiff=iddiff+365
|
|
|
1026 |
enddo
|
|
|
1027 |
iddiff=iddiff+ixday
|
|
|
1028 |
else if (date1(1).lt.date2(1)) then
|
|
|
1029 |
do j=date1(1),date2(1)-1
|
|
|
1030 |
iddiff=iddiff-365
|
|
|
1031 |
enddo
|
|
|
1032 |
iddiff=iddiff-ixday
|
|
|
1033 |
endif
|
|
|
1034 |
|
|
|
1035 |
if (date1(2).gt.date2(2)) then ! compare monthes
|
|
|
1036 |
do j=date2(2),date1(2)-1
|
|
|
1037 |
iddiff=iddiff+idays(j)
|
|
|
1038 |
enddo
|
|
|
1039 |
else if (date1(2).lt.date2(2)) then
|
|
|
1040 |
do j=date1(2),date2(2)-1
|
|
|
1041 |
iddiff=iddiff-idays(j)
|
|
|
1042 |
enddo
|
|
|
1043 |
endif
|
|
|
1044 |
|
|
|
1045 |
iddiff=iddiff+date1(3)-date2(3)
|
|
|
1046 |
ihdiff=iddiff*24+date1(4)-date2(4)
|
|
|
1047 |
imdiff=ihdiff*60+date1(5)-date2(5)
|
|
|
1048 |
|
|
|
1049 |
ihdiff=imdiff/60
|
|
|
1050 |
imdiff=mod(imdiff,60)
|
|
|
1051 |
|
|
|
1052 |
diff=real(ihdiff)+real(imdiff)/60.
|
|
|
1053 |
|
|
|
1054 |
return
|
|
|
1055 |
end
|
|
|
1056 |
subroutine oldtimediff(date1,date2,diff)
|
|
|
1057 |
C ========================================
|
|
|
1058 |
|
|
|
1059 |
C Calculates the time difference in hours for the two dates specified
|
|
|
1060 |
C by the two arrays date1 and date2. They are expected to contain the
|
|
|
1061 |
C following date information:
|
|
|
1062 |
C year month day time step.
|
|
|
1063 |
C
|
|
|
1064 |
C date1 array specifying the first date
|
|
|
1065 |
C date2 array specifying the second date
|
|
|
1066 |
C diff time differenc between date1 and date2 in hours
|
|
|
1067 |
C
|
|
|
1068 |
C Warning: ihdiff is equal to 0 for date1/2 = 880203_12_00 and
|
|
|
1069 |
C 880202_12_24 !!!
|
|
|
1070 |
|
|
|
1071 |
integer date1(5),date2(5)
|
|
|
1072 |
integer idays(12) ! array containing the days of the monthes
|
|
|
1073 |
real diff
|
|
|
1074 |
integer ixday,ihdiff,iddiff,j
|
|
|
1075 |
|
|
|
1076 |
data idays/31,28,31,30,31,30,31,31,30,31,30,31/
|
|
|
1077 |
|
|
|
1078 |
C Determine if the period between date1 and date2 contains a Feb.29
|
|
|
1079 |
|
|
|
1080 |
ixday=0 ! extra day flag
|
|
|
1081 |
|
|
|
1082 |
if (((mod(date1(1),4).eq.0).and.(date1(1).ne.0)).or.
|
|
|
1083 |
> ((mod(date2(1),4).eq.0).and.(date2(1).ne.0))) then
|
|
|
1084 |
if (((date1(2).le.2).and.(date2(2).gt.2)).or.
|
|
|
1085 |
> ((date1(2).gt.2).and.(date2(2).le.2))) then
|
|
|
1086 |
ixday=1 ! set extra day flag to 1
|
|
|
1087 |
idays(2)=29 ! february has 29 days
|
|
|
1088 |
endif
|
|
|
1089 |
endif
|
|
|
1090 |
|
|
|
1091 |
ihdiff=0 ! diff. in hours between date1/date2
|
|
|
1092 |
iddiff=0 ! diff. in days between date1/date2
|
|
|
1093 |
|
|
|
1094 |
if (date1(1).gt.date2(1)) then ! compare years
|
|
|
1095 |
do j=date2(1),date1(1)-1
|
|
|
1096 |
iddiff=iddiff+365+ixday
|
|
|
1097 |
enddo
|
|
|
1098 |
else if (date1(1).lt.date2(1)) then
|
|
|
1099 |
do j=date1(1),date2(1)-1
|
|
|
1100 |
iddiff=iddiff-365-ixday
|
|
|
1101 |
enddo
|
|
|
1102 |
endif
|
|
|
1103 |
|
|
|
1104 |
if (date1(2).gt.date2(2)) then ! compare monthes
|
|
|
1105 |
do j=date2(2),date1(2)-1
|
|
|
1106 |
iddiff=iddiff+idays(j)
|
|
|
1107 |
enddo
|
|
|
1108 |
else if (date1(2).lt.date2(2)) then
|
|
|
1109 |
do j=date1(2),date2(2)-1
|
|
|
1110 |
iddiff=iddiff-idays(j)
|
|
|
1111 |
enddo
|
|
|
1112 |
endif
|
|
|
1113 |
|
|
|
1114 |
iddiff=iddiff+date1(3)-date2(3)
|
|
|
1115 |
ihdiff=iddiff*24+date1(4)-date2(4)+date1(5)-date2(5)
|
|
|
1116 |
|
|
|
1117 |
diff=real(ihdiff)
|
|
|
1118 |
|
|
|
1119 |
return
|
|
|
1120 |
end
|
|
|
1121 |
subroutine newdate(date1,diff,date2)
|
|
|
1122 |
C ====================================
|
|
|
1123 |
C
|
|
|
1124 |
C Routine calculates the new date when diff (in hours) is added to
|
|
|
1125 |
C date1.
|
|
|
1126 |
C date1 int input array contains a date in the form
|
|
|
1127 |
C year,month,day,hour,step
|
|
|
1128 |
C diff real input timestep in hours to go from date1
|
|
|
1129 |
C date2 int output array contains new date in the same form
|
|
|
1130 |
|
|
|
1131 |
integer date1(5),date2(5)
|
|
|
1132 |
integer idays(12) ! array containing the days of the monthes
|
|
|
1133 |
real diff
|
|
|
1134 |
logical yearchange
|
|
|
1135 |
|
|
|
1136 |
data idays/31,28,31,30,31,30,31,31,30,31,30,31/
|
|
|
1137 |
|
|
|
1138 |
yearchange=.false.
|
|
|
1139 |
|
|
|
1140 |
if ((mod(date1(1),4).eq.0).and.(date1(2).le.2)) idays(2)=29
|
|
|
1141 |
|
|
|
1142 |
date2(1)=date1(1)
|
|
|
1143 |
date2(2)=date1(2)
|
|
|
1144 |
date2(3)=date1(3)
|
|
|
1145 |
date2(4)=date1(4)
|
|
|
1146 |
date2(5)=0
|
|
|
1147 |
date2(4)=date1(4)+int(diff)+date1(5)
|
|
|
1148 |
|
|
|
1149 |
if (date2(4).ge.24) then
|
|
|
1150 |
date2(3)=date2(3)+int(date2(4)/24)
|
|
|
1151 |
date2(4)=date2(4)-int(date2(4)/24)*24
|
|
|
1152 |
endif
|
|
|
1153 |
if (date2(4).lt.0) then
|
|
|
1154 |
if (mod(date2(4),24).eq.0) then
|
|
|
1155 |
date2(3)=date2(3)-int(abs(date2(4))/24)
|
|
|
1156 |
date2(4)=date2(4)+int(abs(date2(4))/24)*24
|
|
|
1157 |
else
|
|
|
1158 |
date2(3)=date2(3)-(1+int(abs(date2(4))/24))
|
|
|
1159 |
date2(4)=date2(4)+(1+int(abs(date2(4))/24))*24
|
|
|
1160 |
endif
|
|
|
1161 |
endif
|
|
|
1162 |
|
|
|
1163 |
100 if (date2(3).gt.idays(date2(2))) then
|
|
|
1164 |
if ((date2(2).eq.2).and.(mod(date2(1),4).eq.0)) idays(2)=29
|
|
|
1165 |
date2(3)=date2(3)-idays(date2(2))
|
|
|
1166 |
if (idays(2).eq.29) idays(2)=28
|
|
|
1167 |
date2(2)=date2(2)+1
|
|
|
1168 |
if (date2(2).gt.12) then
|
|
|
1169 |
* date2(1)=date2(1)+int(date2(2)/12)
|
|
|
1170 |
* date2(2)=date2(2)-int(date2(2)/12)*12
|
|
|
1171 |
date2(1)=date2(1)+1
|
|
|
1172 |
date2(2)=date2(2)-12
|
|
|
1173 |
endif
|
|
|
1174 |
if (date2(2).lt.1) then
|
|
|
1175 |
date2(1)=date2(1)-(1+int(abs(date2(2)/12)))
|
|
|
1176 |
date2(2)=date2(2)+(1+int(abs(date2(2)/12)))*12
|
|
|
1177 |
endif
|
|
|
1178 |
goto 100
|
|
|
1179 |
endif
|
|
|
1180 |
200 if (date2(3).lt.1) then
|
|
|
1181 |
date2(2)=date2(2)-1
|
|
|
1182 |
if (date2(2).gt.12) then
|
|
|
1183 |
date2(1)=date2(1)+int(date2(2)/12)
|
|
|
1184 |
date2(2)=date2(2)-int(date2(2)/12)*12
|
|
|
1185 |
endif
|
|
|
1186 |
if (date2(2).lt.1) then
|
|
|
1187 |
date2(1)=date2(1)-(1+int(abs(date2(2)/12)))
|
|
|
1188 |
date2(2)=date2(2)+(1+int(abs(date2(2)/12)))*12
|
|
|
1189 |
endif
|
|
|
1190 |
if ((date2(2).eq.2).and.(mod(date2(1),4).eq.0)) idays(2)=29
|
|
|
1191 |
date2(3)=date2(3)+idays(date2(2))
|
|
|
1192 |
if (idays(2).eq.29) idays(2)=28
|
|
|
1193 |
goto 200
|
|
|
1194 |
endif
|
|
|
1195 |
|
|
|
1196 |
if (date2(2).gt.12) then
|
|
|
1197 |
date2(1)=date2(1)+int(date2(2)/12)
|
|
|
1198 |
date2(2)=date2(2)-int(date2(2)/12)*12
|
|
|
1199 |
endif
|
|
|
1200 |
if (date2(2).lt.1) then
|
|
|
1201 |
date2(1)=date2(1)-(1+int(abs(date2(2)/12)))
|
|
|
1202 |
date2(2)=date2(2)+(1+int(abs(date2(2)/12)))*12
|
|
|
1203 |
endif
|
|
|
1204 |
|
|
|
1205 |
if (date2(1).lt.1000) then
|
|
|
1206 |
if (date2(1).ge.100) date2(1)=date2(1)-100
|
|
|
1207 |
endif
|
|
|
1208 |
|
|
|
1209 |
return
|
|
|
1210 |
end
|
|
|
1211 |
subroutine indipo(lonw,lone,lats,latn,dx,dy,nx,ny,nz,ak,bk,
|
|
|
1212 |
> sp0,sp1,reltpos,levtyp,xpo,ypo,ppo,xind,yind,pind)
|
|
|
1213 |
C =============================================================
|
|
|
1214 |
C
|
|
|
1215 |
C Calculates the position in the grid space from the physical
|
|
|
1216 |
C position.
|
|
|
1217 |
C
|
|
|
1218 |
integer nx,ny,nz,i
|
|
|
1219 |
real dx,dy
|
|
|
1220 |
real ak(nz),bk(nz)
|
|
|
1221 |
real sp0(nx,ny),sp1(nx,ny),spval0,spval1,spval
|
|
|
1222 |
real lonw,lone,lats,latn
|
|
|
1223 |
real xpo,ypo,ppo,reltpos
|
|
|
1224 |
real xind,yind,pind
|
|
|
1225 |
integer levtyp
|
|
|
1226 |
real int2d
|
|
|
1227 |
|
|
|
1228 |
C levtyp=1 pressure levels
|
|
|
1229 |
C levtyp=2 theta levels
|
|
|
1230 |
C levtyp=3 model levels
|
|
|
1231 |
|
|
|
1232 |
C Calculate the grid location to be interpolated to
|
|
|
1233 |
|
|
|
1234 |
xind=(xpo-lonw)/dx+1.
|
|
|
1235 |
yind=(ypo-lats)/dy+1.
|
|
|
1236 |
|
|
|
1237 |
if (nz.eq.1) then
|
|
|
1238 |
pind=1.
|
|
|
1239 |
return
|
|
|
1240 |
endif
|
|
|
1241 |
|
|
|
1242 |
if (levtyp.eq.1) then
|
|
|
1243 |
do i=1,nz-1
|
|
|
1244 |
if ((ak(i).ge.ppo).and.(ak(i+1).le.ppo)) then
|
|
|
1245 |
pind=real(i)+(ak(i)-ppo)/(ak(i)-ak(i+1))
|
|
|
1246 |
return
|
|
|
1247 |
endif
|
|
|
1248 |
enddo
|
|
|
1249 |
else if (levtyp.eq.2) then
|
|
|
1250 |
do i=1,nz-1
|
|
|
1251 |
if ((ak(i).le.ppo).and.(ak(i+1).ge.ppo)) then
|
|
|
1252 |
pind=real(i)+(ak(i)-ppo)/(ak(i)-ak(i+1))
|
|
|
1253 |
return
|
|
|
1254 |
endif
|
|
|
1255 |
enddo
|
|
|
1256 |
else
|
|
|
1257 |
if (ppo.eq.1050.) then
|
|
|
1258 |
pind=1.
|
|
|
1259 |
return
|
|
|
1260 |
else
|
|
|
1261 |
if (reltpos.eq.0.) then
|
|
|
1262 |
spval=int2d(sp0,nx,ny,xind,yind)
|
|
|
1263 |
else if (reltpos.eq.1.) then
|
|
|
1264 |
spval=int2d(sp1,nx,ny,xind,yind)
|
|
|
1265 |
else
|
|
|
1266 |
spval0=int2d(sp0,nx,ny,xind,yind)
|
|
|
1267 |
spval1=int2d(sp1,nx,ny,xind,yind)
|
|
|
1268 |
spval=(1.-reltpos)*spval0+reltpos*spval1
|
|
|
1269 |
endif
|
|
|
1270 |
do i=1,nz-1
|
|
|
1271 |
if ((ak(i)+bk(i)*spval.ge.ppo).and.
|
|
|
1272 |
> (ak(i+1)+bk(i+1)*spval.le.ppo)) then
|
|
|
1273 |
pind=real(i)+(ak(i)+bk(i)*spval-ppo)/
|
|
|
1274 |
> (ak(i)+bk(i)*spval-ak(i+1)-bk(i+1)*spval)
|
|
|
1275 |
return
|
|
|
1276 |
endif
|
|
|
1277 |
enddo
|
|
|
1278 |
endif
|
|
|
1279 |
endif
|
|
|
1280 |
if (levtyp.lt.3) then
|
|
|
1281 |
write(*,*)'error in indipo: ppo ',ppo,' is outside the levels'
|
|
|
1282 |
else
|
|
|
1283 |
if (ppo.gt.(ak(1)+bk(1)*spval)) then
|
|
|
1284 |
* pind=1.
|
|
|
1285 |
pind=(spval-ppo)/(spval-(ak(1)+bk(1)*spval))
|
|
|
1286 |
return
|
|
|
1287 |
else
|
|
|
1288 |
write(*,*)'pos ',xpo,ypo,ppo,' sp ',ak(1)+bk(1)*spval
|
|
|
1289 |
write(*,*)'error in indipo: ppo is outside the levels'
|
|
|
1290 |
endif
|
|
|
1291 |
endif
|
|
|
1292 |
|
|
|
1293 |
return
|
|
|
1294 |
end
|
|
|
1295 |
subroutine nindipo(lonw0,lats0,lonw1,lats1,dx,dy,nx0,ny0,
|
|
|
1296 |
> nx1,ny1,nz,ak,bk,sp0,sp1,reltpos,prelevs,
|
|
|
1297 |
> xpo,ypo,ppo,xind0,yind0,xind1,yind1,pind)
|
|
|
1298 |
C ==============================================================
|
|
|
1299 |
C
|
|
|
1300 |
C Calculates the position in the grid space from the physical
|
|
|
1301 |
C position.
|
|
|
1302 |
C
|
|
|
1303 |
integer nx0,ny0,nx1,ny1,nz,i
|
|
|
1304 |
real dx,dy
|
|
|
1305 |
real ak(nz),bk(nz)
|
|
|
1306 |
real sp0(nx0,ny0),sp1(nx1,ny1),spval0,spval1,spval
|
|
|
1307 |
real lonw0,lats0,lonw1,lats1
|
|
|
1308 |
real xpo,ypo,ppo,reltpos
|
|
|
1309 |
real xind0,yind0,xind1,yind1,pind
|
|
|
1310 |
logical prelevs
|
|
|
1311 |
real int3d
|
|
|
1312 |
|
|
|
1313 |
C Calculate the grid location to be interpolated to
|
|
|
1314 |
|
|
|
1315 |
xind0=(xpo-lonw0)/dx+1.
|
|
|
1316 |
yind0=(ypo-lats0)/dy+1.
|
|
|
1317 |
xind1=(xpo-lonw1)/dx+1.
|
|
|
1318 |
yind1=(ypo-lats1)/dy+1.
|
|
|
1319 |
|
|
|
1320 |
if (prelevs) then
|
|
|
1321 |
do i=1,nz-1
|
|
|
1322 |
if ((ak(i).ge.ppo).and.(ak(i+1).le.ppo)) then
|
|
|
1323 |
pind=real(i)+(ak(i)-ppo)/(ak(i)-ak(i+1))
|
|
|
1324 |
return
|
|
|
1325 |
endif
|
|
|
1326 |
enddo
|
|
|
1327 |
else
|
|
|
1328 |
if (ppo.eq.1050.) then
|
|
|
1329 |
pind=1.
|
|
|
1330 |
return
|
|
|
1331 |
else
|
|
|
1332 |
if (reltpos.eq.0.) then
|
|
|
1333 |
spval=int3d(sp0,nx0,ny0,nz,xind0,yind0,1.)
|
|
|
1334 |
else if (reltpos.eq.1.) then
|
|
|
1335 |
spval=int3d(sp1,nx1,ny1,nz,xind1,yind1,1.)
|
|
|
1336 |
else
|
|
|
1337 |
spval0=int3d(sp0,nx0,ny0,nz,xind0,yind0,1.)
|
|
|
1338 |
spval1=int3d(sp1,nx1,ny1,nz,xind1,yind1,1.)
|
|
|
1339 |
spval=(1.-reltpos)*spval0+reltpos*spval1
|
|
|
1340 |
endif
|
|
|
1341 |
do i=1,nz-1
|
|
|
1342 |
if ((ak(i)+bk(i)*spval.ge.ppo).and.
|
|
|
1343 |
> (ak(i+1)+bk(i+1)*spval.le.ppo)) then
|
|
|
1344 |
pind=real(i)+(ak(i)+bk(i)*spval-ppo)/
|
|
|
1345 |
> (ak(i)+bk(i)*spval-ak(i+1)-bk(i+1)*spval)
|
|
|
1346 |
return
|
|
|
1347 |
endif
|
|
|
1348 |
enddo
|
|
|
1349 |
endif
|
|
|
1350 |
endif
|
|
|
1351 |
if (prelevs) then
|
|
|
1352 |
write(*,*)'error in indipo: ppo is outside the levels'
|
|
|
1353 |
else
|
|
|
1354 |
if (ppo.gt.(ak(1)+bk(1)*spval)) then
|
|
|
1355 |
pind=1.
|
|
|
1356 |
return
|
|
|
1357 |
else
|
|
|
1358 |
write(*,*)'pos ',xpo,ypo,ppo,' sp ',ak(1)+bk(1)*spval
|
|
|
1359 |
write(*,*)'error in indipo: ppo is outside the levels'
|
|
|
1360 |
endif
|
|
|
1361 |
endif
|
|
|
1362 |
|
|
|
1363 |
return
|
|
|
1364 |
end
|
|
|
1365 |
subroutine indipo_mc2(lonw,lone,lats,latn,dx,dy,nx,ny,nz,
|
|
|
1366 |
> bklev,bklay,bkt,zb,
|
|
|
1367 |
> xpo,ypo,ppo,xind,yind,pindlev,pindlay)
|
|
|
1368 |
C =========================================================
|
|
|
1369 |
C
|
|
|
1370 |
C Calculates the position in the grid space from the physical
|
|
|
1371 |
C position.
|
|
|
1372 |
C
|
|
|
1373 |
implicit none
|
|
|
1374 |
|
|
|
1375 |
integer nx,ny,nz,i
|
|
|
1376 |
real dx,dy
|
|
|
1377 |
real bklev(nz),bklay(nz),bkt
|
|
|
1378 |
real zb(nx,ny),spval0,spval1,spval,zbval
|
|
|
1379 |
real lonw,lone,lats,latn
|
|
|
1380 |
real xpo,ypo,ppo
|
|
|
1381 |
real xind,yind,pindlev,pindlay
|
|
|
1382 |
real int3d
|
|
|
1383 |
|
|
|
1384 |
C Calculate the grid location to be interpolated to
|
|
|
1385 |
|
|
|
1386 |
xind=(xpo-lonw)/dx+1.
|
|
|
1387 |
yind=(ypo-lats)/dy+1.
|
|
|
1388 |
|
|
|
1389 |
if (ppo.eq.0.) then
|
|
|
1390 |
pindlev=1.
|
|
|
1391 |
pindlay=1.
|
|
|
1392 |
return
|
|
|
1393 |
endif
|
|
|
1394 |
zbval=int3d(zb,nx,ny,nz,xind,yind,1.)
|
|
|
1395 |
do i=1,nz-1
|
|
|
1396 |
if ((bklev(i)+(1.-bklev(i)/bkt)*zbval.le.ppo).and.
|
|
|
1397 |
> bklev(i+1)+(1.-bklev(i+1)/bkt)*zbval.gt.ppo) then
|
|
|
1398 |
pindlev=real(i)+(ppo-bklev(i)-(1.-bklev(i)/bkt)*zbval)/
|
|
|
1399 |
> ((1.-zbval/bkt)*(bklev(i+1)-bklev(i)))
|
|
|
1400 |
goto 40
|
|
|
1401 |
endif
|
|
|
1402 |
enddo
|
|
|
1403 |
if (ppo.lt.zbval) then
|
|
|
1404 |
write(*,*)'error in indipo_mc2: ppo is below the ground'
|
|
|
1405 |
call exit(1)
|
|
|
1406 |
else if (ppo.lt.bklev(1)+(1.-bklev(1)/bkt)*zbval) then
|
|
|
1407 |
pindlev=(ppo-zbval)/(bklev(1)-bklev(1)/bkt*zbval)
|
|
|
1408 |
return
|
|
|
1409 |
else
|
|
|
1410 |
write(*,*)'pos ',xpo,ypo,ppo,bklev(1)+(1.-bklev(1)/bkt)*zbval
|
|
|
1411 |
write(*,*)'error in indipo_mc2: ppo is outside the levels'
|
|
|
1412 |
call exit(1)
|
|
|
1413 |
endif
|
|
|
1414 |
40 do i=1,nz-1
|
|
|
1415 |
if ((bklay(i)+(1.-bklay(i)/bkt)*zbval.le.ppo).and.
|
|
|
1416 |
> bklay(i+1)+(1.-bklay(i+1)/bkt)*zbval.gt.ppo) then
|
|
|
1417 |
pindlay=real(i)+(ppo-bklay(i)-(1.-bklay(i)/bkt)*zbval)/
|
|
|
1418 |
> ((1.-zbval/bkt)*(bklay(i+1)-bklay(i)))
|
|
|
1419 |
return
|
|
|
1420 |
endif
|
|
|
1421 |
enddo
|
|
|
1422 |
if (ppo.lt.zbval) then
|
|
|
1423 |
write(*,*)'error in indipo_mc2: ppo is below the ground'
|
|
|
1424 |
call exit(1)
|
|
|
1425 |
else if (ppo.lt.bklay(1)+(1.-bklay(1)/bkt)*zbval) then
|
|
|
1426 |
pindlay=(ppo-zbval)/(bklay(1)-bklay(1)/bkt*zbval)
|
|
|
1427 |
return
|
|
|
1428 |
else
|
|
|
1429 |
write(*,*)'pos ',xpo,ypo,ppo,bklay(1)+(1.-bklay(1)/bkt)*zbval
|
|
|
1430 |
write(*,*)'error in indipo_mc2: ppo is outside the layers'
|
|
|
1431 |
call exit(1)
|
|
|
1432 |
endif
|
|
|
1433 |
|
|
|
1434 |
end
|
|
|
1435 |
subroutine ipo(nx,ny,nz,field0,field1,reltpos,
|
|
|
1436 |
> xind,yind,pind,mdv,ipomode,value)
|
|
|
1437 |
C ================================================
|
|
|
1438 |
C
|
|
|
1439 |
C Decides what kind of interpolation should be done.
|
|
|
1440 |
C
|
|
|
1441 |
integer nx,ny,nz
|
|
|
1442 |
real field0(nx,ny,nz)
|
|
|
1443 |
real field1(nx,ny,nz)
|
|
|
1444 |
real reltpos
|
|
|
1445 |
real xind,yind,pind
|
|
|
1446 |
real value,mdv
|
|
|
1447 |
integer ipomode
|
|
|
1448 |
|
|
|
1449 |
if (ipomode.eq.1) then
|
|
|
1450 |
call linipo(nx,ny,nz,field0,field1,reltpos,
|
|
|
1451 |
> xind,yind,pind,mdv,value)
|
|
|
1452 |
else if (ipomode.eq.2) then
|
|
|
1453 |
call bicipo(nx,ny,nz,field0,field1,reltpos,
|
|
|
1454 |
> xind,yind,pind,value)
|
|
|
1455 |
else if (ipomode.eq.3) then
|
|
|
1456 |
call vcuipo(nx,ny,nz,field0,field1,reltpos,
|
|
|
1457 |
> xind,yind,pind,mdv,value)
|
|
|
1458 |
else if (ipomode.eq.11) then
|
|
|
1459 |
call linipolog(nx,ny,nz,field0,field1,reltpos,
|
|
|
1460 |
> xind,yind,pind,mdv,value)
|
|
|
1461 |
else
|
|
|
1462 |
write(*,*)'*** error in SR ipo: unknown ipomode'
|
|
|
1463 |
stop
|
|
|
1464 |
endif
|
|
|
1465 |
|
|
|
1466 |
return
|
|
|
1467 |
end
|
|
|
1468 |
subroutine nipo(nx0,ny0,nx1,ny1,nz,field0,field1,reltpos,
|
|
|
1469 |
> xind0,yind0,xind1,yind1,pind,mdv,ipomode,value)
|
|
|
1470 |
C ==========================================================
|
|
|
1471 |
C
|
|
|
1472 |
C Decides what kind of interpolation should be done.
|
|
|
1473 |
C
|
|
|
1474 |
integer nx0,ny0,nx1,ny1,nz
|
|
|
1475 |
real field0(nx0,ny0,nz)
|
|
|
1476 |
real field1(nx1,ny1,nz)
|
|
|
1477 |
real reltpos
|
|
|
1478 |
real xind0,yind0,xind1,yind1,pind
|
|
|
1479 |
real value,mdv
|
|
|
1480 |
integer ipomode
|
|
|
1481 |
|
|
|
1482 |
if (ipomode.eq.1) then
|
|
|
1483 |
call nlinipo(nx0,ny0,nx1,ny1,nz,field0,field1,reltpos,
|
|
|
1484 |
> xind0,yind0,xind1,yind1,pind,mdv,value)
|
|
|
1485 |
else
|
|
|
1486 |
write(*,*)'*** error in SR ipo: unknown ipomode'
|
|
|
1487 |
stop
|
|
|
1488 |
endif
|
|
|
1489 |
|
|
|
1490 |
return
|
|
|
1491 |
end
|
|
|
1492 |
subroutine linipo(nx,ny,nz,field0,field1,reltpos,
|
|
|
1493 |
> xind,yind,pind,mdv,value)
|
|
|
1494 |
C =================================================
|
|
|
1495 |
C
|
|
|
1496 |
C Does linear interpolation both in time and space.
|
|
|
1497 |
C
|
|
|
1498 |
integer nx,ny,nz
|
|
|
1499 |
real field0(nx,ny,nz)
|
|
|
1500 |
real field1(nx,ny,nz)
|
|
|
1501 |
real reltpos
|
|
|
1502 |
real xind,yind,pind
|
|
|
1503 |
real value,val0,val1,mdv
|
|
|
1504 |
real int3dm,int2dm
|
|
|
1505 |
|
|
|
1506 |
if (reltpos.eq.0.) then
|
|
|
1507 |
if (nz.eq.1) then
|
|
|
1508 |
value=int2dm(field0,nx,ny,xind,yind,mdv)
|
|
|
1509 |
else
|
|
|
1510 |
value=int3dm(field0,nx,ny,nz,xind,yind,pind,mdv)
|
|
|
1511 |
endif
|
|
|
1512 |
else if (reltpos.eq.1.) then
|
|
|
1513 |
if (nz.eq.1) then
|
|
|
1514 |
value=int2dm(field1,nx,ny,xind,yind,mdv)
|
|
|
1515 |
else
|
|
|
1516 |
value=int3dm(field1,nx,ny,nz,xind,yind,pind,mdv)
|
|
|
1517 |
endif
|
|
|
1518 |
else
|
|
|
1519 |
if (nz.eq.1) then
|
|
|
1520 |
val0=int2dm(field0,nx,ny,xind,yind,mdv)
|
|
|
1521 |
val1=int2dm(field1,nx,ny,xind,yind,mdv)
|
|
|
1522 |
else
|
|
|
1523 |
val0=int3dm(field0,nx,ny,nz,xind,yind,pind,mdv)
|
|
|
1524 |
val1=int3dm(field1,nx,ny,nz,xind,yind,pind,mdv)
|
|
|
1525 |
endif
|
|
|
1526 |
if ((val0.eq.mdv).or.(val1.eq.mdv)) then
|
|
|
1527 |
value=mdv
|
|
|
1528 |
else
|
|
|
1529 |
value=(1.-reltpos)*val0+reltpos*val1
|
|
|
1530 |
endif
|
|
|
1531 |
endif
|
|
|
1532 |
|
|
|
1533 |
return
|
|
|
1534 |
end
|
|
|
1535 |
subroutine linipolog(nx,ny,nz,field0,field1,reltpos,
|
|
|
1536 |
> xind,yind,pind,mdv,value)
|
|
|
1537 |
C ====================================================
|
|
|
1538 |
C
|
|
|
1539 |
C Does linear interpolation both in time and space.
|
|
|
1540 |
C
|
|
|
1541 |
integer nx,ny,nz
|
|
|
1542 |
real field0(nx,ny,nz)
|
|
|
1543 |
real field1(nx,ny,nz)
|
|
|
1544 |
real reltpos
|
|
|
1545 |
real xind,yind,pind
|
|
|
1546 |
real value,val0,val1,mdv
|
|
|
1547 |
real int3dmlog
|
|
|
1548 |
|
|
|
1549 |
if (reltpos.eq.0.) then
|
|
|
1550 |
value=int3dmlog(field0,nx,ny,nz,xind,yind,pind,mdv)
|
|
|
1551 |
else if (reltpos.eq.1.) then
|
|
|
1552 |
value=int3dmlog(field1,nx,ny,nz,xind,yind,pind,mdv)
|
|
|
1553 |
else
|
|
|
1554 |
val0=int3dmlog(field0,nx,ny,nz,xind,yind,pind,mdv)
|
|
|
1555 |
val1=int3dmlog(field1,nx,ny,nz,xind,yind,pind,mdv)
|
|
|
1556 |
if ((val0.eq.mdv).or.(val1.eq.mdv)) then
|
|
|
1557 |
value=mdv
|
|
|
1558 |
else
|
|
|
1559 |
value=(1.-reltpos)*val0+reltpos*val1
|
|
|
1560 |
endif
|
|
|
1561 |
endif
|
|
|
1562 |
|
|
|
1563 |
return
|
|
|
1564 |
end
|
|
|
1565 |
subroutine nlinipo(nx0,ny0,nx1,ny1,nz,field0,field1,reltpos,
|
|
|
1566 |
> xind0,yind0,xind1,yind1,pind,mdv,value)
|
|
|
1567 |
C ============================================================
|
|
|
1568 |
C
|
|
|
1569 |
C Does linear interpolation both in time and space.
|
|
|
1570 |
C
|
|
|
1571 |
integer nx0,ny0,nx1,ny1,nz
|
|
|
1572 |
real field0(nx0,ny0,nz)
|
|
|
1573 |
real field1(nx1,ny1,nz)
|
|
|
1574 |
real reltpos
|
|
|
1575 |
real xind0,yind0,xind1,yind1,pind
|
|
|
1576 |
real value,val0,val1,mdv
|
|
|
1577 |
real int3dm
|
|
|
1578 |
|
|
|
1579 |
if (reltpos.eq.0.) then
|
|
|
1580 |
value=int3dm(field0,nx0,ny0,nz,xind0,yind0,pind,mdv)
|
|
|
1581 |
else if (reltpos.eq.1.) then
|
|
|
1582 |
value=int3dm(field1,nx1,ny1,nz,xind1,yind1,pind,mdv)
|
|
|
1583 |
else
|
|
|
1584 |
val0=int3dm(field0,nx0,ny0,nz,xind0,yind0,pind,mdv)
|
|
|
1585 |
val1=int3dm(field1,nx1,ny1,nz,xind1,yind1,pind,mdv)
|
|
|
1586 |
value=(1.-reltpos)*val0+reltpos*val1
|
|
|
1587 |
endif
|
|
|
1588 |
|
|
|
1589 |
return
|
|
|
1590 |
end
|
|
|
1591 |
subroutine vcuipo(nx,ny,nz,field0,field1,reltpos,
|
|
|
1592 |
> xind,yind,pind,mdv,value)
|
|
|
1593 |
C =================================================
|
|
|
1594 |
C
|
|
|
1595 |
C Does linear interpolation both in time and space, except for a cubic
|
|
|
1596 |
C interpolation in the vertical.
|
|
|
1597 |
C
|
|
|
1598 |
integer nx,ny,nz
|
|
|
1599 |
real field0(nx,ny,nz)
|
|
|
1600 |
real field1(nx,ny,nz)
|
|
|
1601 |
real reltpos
|
|
|
1602 |
real xind,yind,pind
|
|
|
1603 |
real value,val0,val1,mdv
|
|
|
1604 |
real int3dmvc
|
|
|
1605 |
|
|
|
1606 |
if (reltpos.eq.0.) then
|
|
|
1607 |
value=int3dmvc(field0,nx,ny,nz,xind,yind,pind,mdv)
|
|
|
1608 |
else if (reltpos.eq.1.) then
|
|
|
1609 |
value=int3dmvc(field1,nx,ny,nz,xind,yind,pind,mdv)
|
|
|
1610 |
else
|
|
|
1611 |
val0=int3dmvc(field0,nx,ny,nz,xind,yind,pind,mdv)
|
|
|
1612 |
val1=int3dmvc(field1,nx,ny,nz,xind,yind,pind,mdv)
|
|
|
1613 |
value=(1.-reltpos)*val0+reltpos*val1
|
|
|
1614 |
endif
|
|
|
1615 |
|
|
|
1616 |
return
|
|
|
1617 |
end
|
|
|
1618 |
subroutine linipom(nx,ny,nz,field0,field1,reltpos,
|
|
|
1619 |
> mdv,xind,yind,pind,value)
|
|
|
1620 |
C =================================================
|
|
|
1621 |
C
|
|
|
1622 |
C Does linear interpolation both in time and space.
|
|
|
1623 |
C
|
|
|
1624 |
integer nx,ny,nz
|
|
|
1625 |
real field0(nx,ny,nz)
|
|
|
1626 |
real field1(nx,ny,nz)
|
|
|
1627 |
real reltpos,mdv
|
|
|
1628 |
real xind,yind,pind
|
|
|
1629 |
real value,val0,val1
|
|
|
1630 |
real int3dm
|
|
|
1631 |
|
|
|
1632 |
if (reltpos.eq.0.) then
|
|
|
1633 |
value=int3dm(field0,nx,ny,nz,xind,yind,pind,mdv)
|
|
|
1634 |
else if (reltpos.eq.1.) then
|
|
|
1635 |
value=int3dm(field1,nx,ny,nz,xind,yind,pind,mdv)
|
|
|
1636 |
else
|
|
|
1637 |
val0=int3dm(field0,nx,ny,nz,xind,yind,pind,mdv)
|
|
|
1638 |
val1=int3dm(field1,nx,ny,nz,xind,yind,pind,mdv)
|
|
|
1639 |
value=(1.-reltpos)*val0+reltpos*val1
|
|
|
1640 |
endif
|
|
|
1641 |
|
|
|
1642 |
return
|
|
|
1643 |
end
|
|
|
1644 |
subroutine bicipo(nx,ny,nz,fld0,fld1,reltpos,
|
|
|
1645 |
> xind,yind,pind,value)
|
|
|
1646 |
C ===================================================
|
|
|
1647 |
C
|
|
|
1648 |
C Does linear interpolation in time and bicubic interpolation
|
|
|
1649 |
C in space (horizontally)
|
|
|
1650 |
C
|
|
|
1651 |
integer nx,ny,nz
|
|
|
1652 |
real fld0(nx,ny,nz),fld1(nx,ny,nz)
|
|
|
1653 |
real reltpos,value
|
|
|
1654 |
real xind,yind,pind
|
|
|
1655 |
|
|
|
1656 |
C Declarations of local variables
|
|
|
1657 |
|
|
|
1658 |
integer k1,k2
|
|
|
1659 |
real frac1k,frac2k
|
|
|
1660 |
real val0,val1,val10,val20,val11,val21
|
|
|
1661 |
real cint2d
|
|
|
1662 |
|
|
|
1663 |
k1=int(pind)
|
|
|
1664 |
k2=k1+1
|
|
|
1665 |
|
|
|
1666 |
frac1k=pind-aint(pind)
|
|
|
1667 |
frac2k=1.-frac1k
|
|
|
1668 |
|
|
|
1669 |
C Interpolate value on lower level for first time
|
|
|
1670 |
|
|
|
1671 |
if ((reltpos.ne.1.).and.(frac2k.ne.0.)) then
|
|
|
1672 |
val10=cint2d(fld0,nx,ny,nz,xind,yind,k1)
|
|
|
1673 |
endif
|
|
|
1674 |
|
|
|
1675 |
C Interpolate value on upper level for first time
|
|
|
1676 |
|
|
|
1677 |
if ((reltpos.ne.1.).and.(frac1k.ne.0.)) then
|
|
|
1678 |
val20=cint2d(fld0,nx,ny,nz,xind,yind,k2)
|
|
|
1679 |
endif
|
|
|
1680 |
|
|
|
1681 |
C Interpolate value on lower level for second time
|
|
|
1682 |
|
|
|
1683 |
if ((reltpos.ne.0.).and.(frac2k.ne.0.)) then
|
|
|
1684 |
val11=cint2d(fld1,nx,ny,nz,xind,yind,k1)
|
|
|
1685 |
endif
|
|
|
1686 |
|
|
|
1687 |
C Interpolate value on upper level for second time
|
|
|
1688 |
|
|
|
1689 |
if ((reltpos.ne.0.).and.(frac1k.ne.0.)) then
|
|
|
1690 |
val21=cint2d(fld1,nx,ny,nz,xind,yind,k2)
|
|
|
1691 |
endif
|
|
|
1692 |
|
|
|
1693 |
C Do vertical interpolation and interpolation in time
|
|
|
1694 |
|
|
|
1695 |
if (reltpos.eq.0.) then
|
|
|
1696 |
value=frac2k*val10+frac1k*val20
|
|
|
1697 |
else if (reltpos.eq.1.) then
|
|
|
1698 |
value=frac2k*val11+frac1k*val21
|
|
|
1699 |
else
|
|
|
1700 |
val0=frac2k*val10+frac1k*val20
|
|
|
1701 |
val1=frac2k*val11+frac1k*val21
|
|
|
1702 |
value=(1.-reltpos)*val0+reltpos*val1
|
|
|
1703 |
endif
|
|
|
1704 |
|
|
|
1705 |
return
|
|
|
1706 |
end
|