| 3 |
michaesp |
1 |
PROGRAM select
|
|
|
2 |
|
|
|
3 |
c **************************************************************
|
|
|
4 |
c * Select trajectories from LSL file *
|
|
|
5 |
c * Michael Sprenger / January, February 2008 *
|
|
|
6 |
c **************************************************************
|
|
|
7 |
|
|
|
8 |
implicit none
|
|
|
9 |
|
|
|
10 |
c --------------------------------------------------------------
|
|
|
11 |
c Declaration of parameters
|
|
|
12 |
c --------------------------------------------------------------
|
|
|
13 |
|
|
|
14 |
c Maximum number of columns per trajectory
|
|
|
15 |
integer maxcol
|
|
|
16 |
parameter (maxcol=100)
|
|
|
17 |
|
|
|
18 |
c Maximum number of commands
|
|
|
19 |
integer maxcmd
|
|
|
20 |
parameter (maxcmd=10000)
|
|
|
21 |
|
|
|
22 |
c --------------------------------------------------------------
|
|
|
23 |
c Declaration of variables
|
|
|
24 |
c --------------------------------------------------------------
|
|
|
25 |
|
|
|
26 |
c Input and output files
|
|
|
27 |
character*120 inp_lslfile ! Input lsl file
|
|
|
28 |
character*120 out_lslfile ! Output lsl file
|
|
|
29 |
character*120 inp_criteria ! Input file with criteria
|
|
|
30 |
integer inpmode ! Format of input file
|
|
|
31 |
integer outmode ! Format of output file
|
|
|
32 |
character*80 outformat ! Trajectory/Boolean/Index of output
|
|
|
33 |
character*80 regionf ! Name of the regionfile
|
|
|
34 |
|
|
|
35 |
c Trajectories
|
|
|
36 |
integer ntim ! Number of times
|
|
|
37 |
integer ncol ! Number of columns (including time...)
|
|
|
38 |
integer ntrainp,ntraout ! Number of trajectories
|
|
|
39 |
character*80 vars(maxcol) ! Names of trajectory columns
|
|
|
40 |
integer basetime(6) ! Base time of trajectory (first line in lsl)
|
|
|
41 |
real,allocatable, dimension (:,:,:) :: trainp ! Input trajectories
|
|
|
42 |
real,allocatable, dimension (:,:,:) :: traout ! Output trajectories
|
|
|
43 |
integer,allocatable,dimension (:) :: selflag ! Flag for selection
|
|
|
44 |
real,allocatable, dimension (:) :: time ! Times of the trajectory
|
|
|
45 |
integer,allocatable,dimension (:,:) :: trigger ! Trigger column
|
|
|
46 |
integer itrigger ! Column index for trigger
|
|
|
47 |
character*80 addtrigger ! Flag whether to add trigger column
|
|
|
48 |
|
|
|
49 |
c Command stack
|
|
|
50 |
real cmd(maxcmd) ! Decoded slection criterion
|
|
|
51 |
integer ncmd ! Number of commands
|
|
|
52 |
|
|
|
53 |
c Common block for initialisation of polygon check
|
|
|
54 |
real tlonv(2000),vlat_c(2000),vlon_c(2000)
|
|
|
55 |
real xlat_c,xlon_c
|
|
|
56 |
integer ibndry,nv_c
|
|
|
57 |
data ibndry/0/
|
|
|
58 |
common /spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry
|
|
|
59 |
|
|
|
60 |
c Auxiliary variables
|
|
|
61 |
integer stat ! Logical (state) variable
|
|
|
62 |
integer fid,fod,fcr ! File identifier for input and output
|
|
|
63 |
integer i,j,k ! Index counter
|
|
|
64 |
integer isok ! Flag for selected trajectory
|
|
|
65 |
real param(1000) ! List of parameters
|
|
|
66 |
integer nparam ! Number of parameters
|
|
|
67 |
character*80 specialstr ! Name of special command
|
|
|
68 |
integer len
|
|
|
69 |
integer,allocatable,dimension (:) :: trigger1 ! Trigger column
|
|
|
70 |
real,allocatable,dimension (:,:) :: trainp1 ! A single trajectory
|
|
|
71 |
character ch
|
|
|
72 |
|
|
|
73 |
c --------------------------------------------------------------
|
|
|
74 |
c Preparations
|
|
|
75 |
c --------------------------------------------------------------
|
|
|
76 |
|
|
|
77 |
c Write start message
|
|
|
78 |
print*,'========================================================='
|
|
|
79 |
print*,' *** START OF PROGRAM SELECT ***'
|
|
|
80 |
print*
|
|
|
81 |
|
|
|
82 |
c Read parameter file
|
|
|
83 |
open(10,file='select.param')
|
|
|
84 |
read(10,*) inp_lslfile
|
|
|
85 |
read(10,*) out_lslfile
|
|
|
86 |
read(10,*) outformat
|
|
|
87 |
read(10,*) inp_criteria
|
|
|
88 |
read(10,*) ntrainp
|
|
|
89 |
read(10,*) ntim
|
|
|
90 |
read(10,*) ncol
|
|
|
91 |
read(10,*) regionf
|
|
|
92 |
read(10,*) addtrigger
|
|
|
93 |
close(10)
|
|
|
94 |
|
|
|
95 |
c Set the formats of the input and output files
|
|
|
96 |
call mode_tra(inpmode,inp_lslfile)
|
|
|
97 |
if (inpmode.eq.-1) inpmode=1
|
|
|
98 |
call mode_tra(outmode,out_lslfile)
|
|
|
99 |
if ( (outmode.eq.-1).and.(outformat.ne.'startf') ) then
|
|
|
100 |
outmode=1
|
|
|
101 |
endif
|
|
|
102 |
|
|
|
103 |
c Allocate memory for a single trajectory
|
|
|
104 |
allocate(trainp(ntrainp,ntim,ncol),stat=stat)
|
|
|
105 |
if (stat.ne.0) stop '*** error allocating array trainp ***'
|
|
|
106 |
allocate(time(ntim),stat=stat)
|
|
|
107 |
if (stat.ne.0) stop '*** error allocating array time ***'
|
|
|
108 |
allocate(selflag(ntrainp),stat=stat)
|
|
|
109 |
if (stat.ne.0) stop '*** error allocating array selflag ***'
|
|
|
110 |
allocate(trigger(ntrainp,ntim),stat=stat)
|
|
|
111 |
if (stat.ne.0) stop '*** error allocating array trigger ***'
|
|
|
112 |
allocate(trigger1(ntim),stat=stat)
|
|
|
113 |
if (stat.ne.0) stop '*** error allocating array trigger1 ***'
|
|
|
114 |
allocate(trainp1(ntim,ncol),stat=stat)
|
|
|
115 |
if (stat.ne.0) stop '*** error allocating array trainp1 ***'
|
|
|
116 |
|
|
|
117 |
c Read the input trajectory file
|
|
|
118 |
call ropen_tra(fid,inp_lslfile,ntrainp,ntim,ncol,
|
|
|
119 |
> basetime,vars,inpmode)
|
|
|
120 |
call read_tra (fid,trainp,ntrainp,ntim,ncol,inpmode)
|
|
|
121 |
call close_tra(fid,inpmode)
|
|
|
122 |
|
|
|
123 |
c Check that first four columns correspond to time,lon,lat,p
|
|
|
124 |
if ( (vars(1).ne.'time' ).or.
|
|
|
125 |
> (vars(2).ne.'xpos' ).and.(vars(2).ne.'lon' ).or.
|
|
|
126 |
> (vars(3).ne.'ypos' ).and.(vars(3).ne.'lat' ).or.
|
|
|
127 |
> (vars(4).ne.'ppos' ).and.(vars(4).ne.'p' ) )
|
|
|
128 |
>then
|
|
|
129 |
print*,' ERROR: problem with input trajectories ...'
|
|
|
130 |
stop
|
|
|
131 |
endif
|
|
|
132 |
vars(1) = 'time'
|
|
|
133 |
vars(2) = 'lon'
|
|
|
134 |
vars(3) = 'lat'
|
|
|
135 |
vars(4) = 'p'
|
|
|
136 |
|
|
|
137 |
c Get the trajectory times from first trajectory
|
|
|
138 |
do i=1,ntim
|
|
|
139 |
time(i)=trainp(1,i,1)
|
|
|
140 |
enddo
|
|
|
141 |
|
|
|
142 |
c Init the trigger field - first check whether it is already available
|
|
|
143 |
itrigger = 0
|
|
|
144 |
do i=1,ncol
|
|
|
145 |
if ( vars(i).eq.'TRIGGER' ) itrigger = i
|
|
|
146 |
enddo
|
|
|
147 |
|
|
|
148 |
if ( itrigger.ne.0 ) then
|
|
|
149 |
do i=1,ntrainp
|
|
|
150 |
do j=1,ntim
|
|
|
151 |
trigger(i,j) = nint( trainp(i,j,itrigger) )
|
|
|
152 |
enddo
|
|
|
153 |
enddo
|
|
|
154 |
else
|
|
|
155 |
do i=1,ntrainp
|
|
|
156 |
do j=1,ntim
|
|
|
157 |
trigger(i,j) = 0
|
|
|
158 |
enddo
|
|
|
159 |
enddo
|
|
|
160 |
endif
|
|
|
161 |
|
|
|
162 |
c Write some info about the trajectory
|
|
|
163 |
print*,'---- INPUT PARAMETERS -----------------------------------'
|
|
|
164 |
write(*,*)
|
|
|
165 |
write(*,*) 'Input file : ',trim(inp_lslfile)
|
|
|
166 |
write(*,*) 'Output file : ',trim(out_lslfile)
|
|
|
167 |
write(*,*) 'Output format : ',trim(outformat)
|
|
|
168 |
write(*,*) 'Criteria file : ',trim(inp_criteria)
|
|
|
169 |
write(*,*) '# tra : ',ntrainp
|
|
|
170 |
write(*,*) '# time : ',ntim
|
|
|
171 |
write(*,*) '# col : ',ncol
|
|
|
172 |
write(*,*) 'Region file : ',trim(regionf)
|
|
|
173 |
write(*,*) 'Add trigger : ',trim(addtrigger)
|
|
|
174 |
|
|
|
175 |
print*
|
|
|
176 |
print*,'---- INPUT TRAJECTORY FILE ------------------------------'
|
|
|
177 |
print*
|
|
|
178 |
write(*,'(1x,a12,i4,a10)') 'Vars : ',1,trim(vars(1))
|
|
|
179 |
do i=2,ncol
|
|
|
180 |
write(*,'(1x,a12,i4,a10)') ' ',i,trim(vars(i))
|
|
|
181 |
enddo
|
|
|
182 |
print*
|
|
|
183 |
write(*,'(1x,a12,i4,f10.2)') 'Time : ',1,time(1)
|
|
|
184 |
do i=2,ntim
|
|
|
185 |
write(*,'(1x,a12,i4,f10.2)') ' ',i,time(i)
|
|
|
186 |
enddo
|
|
|
187 |
print*
|
|
|
188 |
write(*,'(1x,a12,i4,i10)') 'Base date : ',1,basetime(1)
|
|
|
189 |
write(*,'(1x,a12,i4,i10)') ' ',2,basetime(2)
|
|
|
190 |
write(*,'(1x,a12,i4,i10)') ' ',3,basetime(3)
|
|
|
191 |
write(*,'(1x,a12,i4,i10)') ' ',4,basetime(4)
|
|
|
192 |
write(*,'(1x,a12,i4,i10)') ' ',5,basetime(5)
|
|
|
193 |
print*
|
|
|
194 |
write(*,'(1x,a12,i4,i10)') 'Time range : ',6,basetime(6)
|
|
|
195 |
print*
|
|
|
196 |
if ( itrigger.ne.0 ) then
|
|
|
197 |
print*,'TRIGGER FIELD FOUND IN COLUMN ',itrigger
|
|
|
198 |
print*
|
|
|
199 |
endif
|
|
|
200 |
|
|
|
201 |
c Read and decode the selection criterion
|
|
|
202 |
fcr = 10
|
|
|
203 |
open(fcr,file=inp_criteria)
|
|
|
204 |
ncmd=maxcmd
|
|
|
205 |
call decode(fcr,cmd,ncmd,vars,ncol,time,ntim,regionf)
|
|
|
206 |
close(fcr)
|
|
|
207 |
|
|
|
208 |
print*
|
|
|
209 |
print*,'---- PSEUDO CODE FOR SELECTION --------------------------'
|
|
|
210 |
print*
|
|
|
211 |
call dumpcode(cmd,ncmd,vars,ncol,time,ntim)
|
|
|
212 |
|
|
|
213 |
c --------------------------------------------------------------
|
|
|
214 |
c Loop over all trajectories - selection
|
|
|
215 |
c --------------------------------------------------------------
|
|
|
216 |
|
|
|
217 |
c Prepare string and parameters for SPECIAL commands
|
|
|
218 |
if ( cmd(1).eq.0 ) then
|
|
|
219 |
|
|
|
220 |
c Get command string
|
|
|
221 |
j = 2
|
|
|
222 |
len = nint(cmd(j))
|
|
|
223 |
specialstr = ''
|
|
|
224 |
do k=1,len
|
|
|
225 |
j = j + 1
|
|
|
226 |
specialstr = trim(specialstr)//char(nint(cmd(j)))
|
|
|
227 |
enddo
|
|
|
228 |
|
|
|
229 |
c Get paramters
|
|
|
230 |
j = j + 1
|
|
|
231 |
nparam = nint(cmd(j))
|
|
|
232 |
do k=1,nparam
|
|
|
233 |
j = j + 1
|
|
|
234 |
param(k) = cmd(j)
|
|
|
235 |
enddo
|
|
|
236 |
|
|
|
237 |
endif
|
|
|
238 |
|
|
|
239 |
c Init the counter for selected trajectories
|
|
|
240 |
ntraout = 0
|
|
|
241 |
|
|
|
242 |
c Loop over all trajectories
|
|
|
243 |
do i=1,ntrainp
|
|
|
244 |
|
|
|
245 |
c Decide whether the trajectory is selected (handle SPECIAL commands)
|
|
|
246 |
if (cmd(1).ne.0 ) then
|
|
|
247 |
|
|
|
248 |
c Copy a single trajectory to <trainp1> and <trigger1>
|
|
|
249 |
do j=1,ntim
|
|
|
250 |
do k=1,ncol
|
|
|
251 |
trainp1(j,k) = trainp(i,j,k)
|
|
|
252 |
enddo
|
|
|
253 |
trigger1(j) = trigger(i,j)
|
|
|
254 |
enddo
|
|
|
255 |
|
|
|
256 |
c Do the selection
|
|
|
257 |
call select_tra (isok,cmd,ncmd,trainp1,trigger1,ntim,ncol)
|
|
|
258 |
|
|
|
259 |
c The trigger might be changed in the selection - copy it
|
|
|
260 |
do j=1,ntim
|
|
|
261 |
trigger(i,j) = trigger1(j)
|
|
|
262 |
enddo
|
|
|
263 |
|
|
|
264 |
else
|
|
|
265 |
call special (isok,specialstr,trainp(i,:,:),ntim,ncol,
|
|
|
266 |
> vars,time,param,nparam)
|
|
|
267 |
endif
|
|
|
268 |
|
|
|
269 |
c Set flag for selected trajectories
|
|
|
270 |
if (isok.eq.1) then
|
|
|
271 |
selflag(i) = 1
|
|
|
272 |
ntraout = ntraout + 1
|
|
|
273 |
else
|
|
|
274 |
selflag(i) = 0
|
|
|
275 |
endif
|
|
|
276 |
|
|
|
277 |
enddo
|
|
|
278 |
|
|
|
279 |
c --------------------------------------------------------------
|
|
|
280 |
c Write output trajectories
|
|
|
281 |
c --------------------------------------------------------------
|
|
|
282 |
|
|
|
283 |
c ------ Write output trajectories -----------------------------
|
|
|
284 |
if ( outformat.eq.'trajectory' ) then
|
|
|
285 |
|
|
|
286 |
|
|
|
287 |
c Define the trigger field if it is not yet defined
|
|
|
288 |
if ( ( addtrigger.eq.'-trigger' ).and.(itrigger.eq.0) ) then
|
|
|
289 |
ncol = ncol + 1
|
|
|
290 |
vars(ncol) = 'TRIGGER'
|
|
|
291 |
itrigger = ncol
|
|
|
292 |
endif
|
|
|
293 |
|
|
|
294 |
c Allocate memory for output trajectory
|
|
|
295 |
allocate(traout(ntraout,ntim,ncol),stat=stat)
|
|
|
296 |
if (stat.ne.0) stop '*** error allocating array apply ***'
|
|
|
297 |
|
|
|
298 |
c Set output trajectories
|
|
|
299 |
j = 0
|
|
|
300 |
do i=1,ntrainp
|
|
|
301 |
if (selflag(i).eq.1) then
|
|
|
302 |
j = j + 1
|
|
|
303 |
traout(j,1:ntim,1:ncol) = trainp(i,1:ntim,1:ncol)
|
|
|
304 |
if ( itrigger.ne.0 ) then
|
|
|
305 |
traout(j,1:ntim,itrigger) = real(trigger(i,1:ntim))
|
|
|
306 |
endif
|
|
|
307 |
endif
|
|
|
308 |
enddo
|
|
|
309 |
|
|
|
310 |
c Write trajectories
|
|
|
311 |
call wopen_tra(fod,out_lslfile,ntraout,ntim,ncol,
|
|
|
312 |
> basetime,vars,outmode)
|
|
|
313 |
call write_tra(fod,traout,ntraout,ntim,ncol,outmode)
|
|
|
314 |
call close_tra(fod,outmode)
|
|
|
315 |
|
|
|
316 |
c ------ Write index list -------------------------------------
|
|
|
317 |
elseif ( outformat.eq.'index' ) then
|
|
|
318 |
|
|
|
319 |
open(10,file=out_lslfile)
|
|
|
320 |
do i=1,ntrainp
|
|
|
321 |
if ( selflag(i).eq.1) write(10,*) i
|
|
|
322 |
enddo
|
|
|
323 |
close(10)
|
|
|
324 |
|
|
|
325 |
c ------ Write boolean list -----------------------------------
|
|
|
326 |
elseif ( outformat.eq.'boolean' ) then
|
|
|
327 |
|
|
|
328 |
open(10,file=out_lslfile)
|
|
|
329 |
do i=1,ntrainp
|
|
|
330 |
write(10,'(i1)') selflag(i)
|
|
|
331 |
enddo
|
|
|
332 |
close(10)
|
|
|
333 |
|
|
|
334 |
c ------ Write count -------------------------------------------
|
|
|
335 |
elseif ( outformat.eq.'count' ) then
|
|
|
336 |
|
|
|
337 |
open(10,file=out_lslfile)
|
|
|
338 |
write(10,'(i7)') ntraout
|
|
|
339 |
close(10)
|
|
|
340 |
|
|
|
341 |
c ------ Write starting positions -----------------------------
|
|
|
342 |
elseif ( outformat.eq.'startf' ) then
|
|
|
343 |
|
|
|
344 |
c Allocate memory for output trajectory
|
|
|
345 |
allocate(traout(ntraout,1,ncol),stat=stat)
|
|
|
346 |
if (stat.ne.0) stop '*** error allocating array apply ***'
|
|
|
347 |
|
|
|
348 |
c Set output trajectories
|
|
|
349 |
j = 0
|
|
|
350 |
do i=1,ntrainp
|
|
|
351 |
if (selflag(i).eq.1) then
|
|
|
352 |
j = j + 1
|
|
|
353 |
traout(j,1,:) = trainp(i,1,:)
|
|
|
354 |
endif
|
|
|
355 |
enddo
|
|
|
356 |
|
|
|
357 |
c Write trajectories
|
|
|
358 |
if (outmode.ne.-1) then
|
|
|
359 |
|
|
|
360 |
call wopen_tra(fod,out_lslfile,ntraout,1,ncol,
|
|
|
361 |
> basetime,vars,outmode)
|
|
|
362 |
call write_tra(fod,traout,ntraout,1,ncol,outmode)
|
|
|
363 |
call close_tra(fod,outmode)
|
|
|
364 |
|
|
|
365 |
c Output as a triple list (corresponding to <startf> file)
|
|
|
366 |
else
|
|
|
367 |
|
|
|
368 |
fid = 10
|
|
|
369 |
open(fid,file=out_lslfile)
|
|
|
370 |
do i=1,ntraout
|
|
|
371 |
write(fid,'(3f10.3)') traout(i,1,2), ! longitude
|
|
|
372 |
> traout(i,1,3), ! latitude
|
|
|
373 |
> traout(i,1,4) ! pressure
|
|
|
374 |
enddo
|
|
|
375 |
close(fid)
|
|
|
376 |
|
|
|
377 |
endif
|
|
|
378 |
|
|
|
379 |
endif
|
|
|
380 |
|
|
|
381 |
c Write some status information, and end of program message
|
|
|
382 |
print*
|
|
|
383 |
print*,'---- STATUS INFORMATION --------------------------------'
|
|
|
384 |
print*
|
|
|
385 |
print*,' # input trajectories : ',ntrainp
|
|
|
386 |
print*,' # output trajectories : ',ntraout
|
|
|
387 |
print*
|
|
|
388 |
print*,' *** END OF PROGRAM SELECT ***'
|
|
|
389 |
print*,'========================================================='
|
|
|
390 |
|
|
|
391 |
stop
|
|
|
392 |
|
|
|
393 |
c Exception handling
|
|
|
394 |
100 stop 'select: First column in input trajectory must be <time>'
|
|
|
395 |
101 stop 'select: Input trajectory file is empty'
|
|
|
396 |
|
|
|
397 |
end
|
|
|
398 |
|
|
|
399 |
c --------------------------------------------------------------
|
|
|
400 |
c Dump the command list
|
|
|
401 |
c --------------------------------------------------------------
|
|
|
402 |
|
|
|
403 |
subroutine dumpcode(out,n,vars,nvars,times,ntimes)
|
|
|
404 |
|
|
|
405 |
c Write the command list to screen. The command list is decoded
|
|
|
406 |
c by call to <decode>
|
|
|
407 |
|
|
|
408 |
implicit none
|
|
|
409 |
|
|
|
410 |
c Declaration of subroutine parameters
|
|
|
411 |
integer n
|
|
|
412 |
real out(n)
|
|
|
413 |
integer nvars
|
|
|
414 |
character*80 vars(nvars)
|
|
|
415 |
integer ntimes
|
|
|
416 |
real times(ntimes)
|
|
|
417 |
|
|
|
418 |
c A single command
|
|
|
419 |
character*80 cmd
|
|
|
420 |
character*80 var,mode,strtim
|
|
|
421 |
integer nval
|
|
|
422 |
integer ntim
|
|
|
423 |
integer ivar,imode,icmd,itime
|
|
|
424 |
|
|
|
425 |
c Auxiliary variables
|
|
|
426 |
integer i,j
|
|
|
427 |
|
|
|
428 |
c Loop through the complete list
|
|
|
429 |
i=0
|
|
|
430 |
100 if (i.lt.n) then
|
|
|
431 |
|
|
|
432 |
write(*,*) '---------------------------------------'
|
|
|
433 |
|
|
|
434 |
c Get command
|
|
|
435 |
i=i+1
|
|
|
436 |
icmd=nint(out(i))
|
|
|
437 |
|
|
|
438 |
c Special handling of SPECIAL commands
|
|
|
439 |
if ( icmd.eq.0 ) then
|
|
|
440 |
|
|
|
441 |
c Write 'header' for SPECIAL command
|
|
|
442 |
write(*,'(i5,f15.4,10x,a10)') i,out(i),'SPECIAL'
|
|
|
443 |
|
|
|
444 |
c Write command string
|
|
|
445 |
i = i + 1
|
|
|
446 |
icmd = nint(out(i))
|
|
|
447 |
write(*,'(i5,f15.4,10x,a10)') i,out(i),'LEN(CMD)'
|
|
|
448 |
do j=1,icmd
|
|
|
449 |
i = i + 1
|
|
|
450 |
ivar = nint(out(i))
|
|
|
451 |
write(*,'(i5,f15.4,10x,a10)') i,out(i),char(ivar)
|
|
|
452 |
enddo
|
|
|
453 |
|
|
|
454 |
c Write parameters
|
|
|
455 |
i = i + 1
|
|
|
456 |
nval = nint(out(i))
|
|
|
457 |
write(*,'(i5,f15.4,10x,a10)') i,out(i),'#PARAMETER'
|
|
|
458 |
do j=1,nval
|
|
|
459 |
i=i+1
|
|
|
460 |
if ( var.ne.'INPOLYGON' ) then
|
|
|
461 |
write(*,'(i5,f15.4)') i,out(i)
|
|
|
462 |
else
|
|
|
463 |
write(*,'(i5,f15.4,a2)') i,out(i),char(nint(out(i)))
|
|
|
464 |
endif
|
|
|
465 |
enddo
|
|
|
466 |
|
|
|
467 |
c Nothing else to do - exit
|
|
|
468 |
goto 200
|
|
|
469 |
|
|
|
470 |
endif
|
|
|
471 |
|
|
|
472 |
c Set the command
|
|
|
473 |
if (icmd.eq. 1) cmd='GT'
|
|
|
474 |
if (icmd.eq. 2) cmd='LT'
|
|
|
475 |
if (icmd.eq. 3) cmd='IN'
|
|
|
476 |
if (icmd.eq. 4) cmd='OUT'
|
|
|
477 |
if (icmd.eq. 5) cmd='EQ'
|
|
|
478 |
if (icmd.eq. 6) cmd='TRUE'
|
|
|
479 |
if (icmd.eq. 7) cmd='FALSE'
|
|
|
480 |
if (icmd.eq. 8) cmd='ALL'
|
|
|
481 |
if (icmd.eq. 9) cmd='ANY'
|
|
|
482 |
if (icmd.eq. 10) cmd='NONE'
|
|
|
483 |
if (icmd.eq. -1) cmd='BEGIN'
|
|
|
484 |
if (icmd.eq. -2) cmd='END'
|
|
|
485 |
if (icmd.eq. -3) cmd='AND'
|
|
|
486 |
if (icmd.eq. -4) cmd='OR'
|
|
|
487 |
write(*,'(i5,f15.4,10x,a10)') i,out(i),trim(cmd)
|
|
|
488 |
if (icmd.lt.0) goto 100
|
|
|
489 |
|
|
|
490 |
c Get variable
|
|
|
491 |
i=i+1
|
|
|
492 |
ivar=nint(out(i))
|
|
|
493 |
if ( ivar.eq. -1 ) then
|
|
|
494 |
var = 'DIST'
|
|
|
495 |
elseif ( ivar.eq. -2 ) then
|
|
|
496 |
var = 'DIST0'
|
|
|
497 |
elseif ( ivar.eq. -3 ) then
|
|
|
498 |
var = 'INPOLYGON'
|
|
|
499 |
elseif ( ivar.eq. -4 ) then
|
|
|
500 |
var = 'INBOX'
|
|
|
501 |
elseif ( ivar.eq. -5 ) then
|
|
|
502 |
var = 'INCIRCLE'
|
|
|
503 |
elseif ( ivar.eq. -6 ) then
|
|
|
504 |
var = 'INREGION'
|
|
|
505 |
elseif ( ivar.eq. -7 ) then
|
|
|
506 |
var = 'TRIGGER'
|
|
|
507 |
else
|
|
|
508 |
var=vars(ivar)
|
|
|
509 |
endif
|
|
|
510 |
write(*,'(i5,f15.4,10x,a10)') i,out(i),trim(var)
|
|
|
511 |
|
|
|
512 |
c Get variable mode
|
|
|
513 |
i=i+1
|
|
|
514 |
imode=nint(out(i))
|
|
|
515 |
if (imode.eq.1) mode='VALUE'
|
|
|
516 |
if (imode.eq.2) mode='MEAN'
|
|
|
517 |
if (imode.eq.3) mode='MAX'
|
|
|
518 |
if (imode.eq.4) mode='MIN'
|
|
|
519 |
if (imode.eq.5) mode='VAR'
|
|
|
520 |
if (imode.eq.6) mode='SUM'
|
|
|
521 |
if (imode.eq.7) mode='CHANGE'
|
|
|
522 |
if (imode.eq.8) mode='DIFF'
|
|
|
523 |
write(*,'(i5,f15.4,10x,a10)') i,out(i),trim(mode)
|
|
|
524 |
|
|
|
525 |
c Get values
|
|
|
526 |
i=i+1
|
|
|
527 |
nval=nint(out(i))
|
|
|
528 |
write(*,'(i5,f15.4,10x,a10)') i,out(i),'#PARAMETER'
|
|
|
529 |
do j=1,nval
|
|
|
530 |
i=i+1
|
|
|
531 |
|
|
|
532 |
if ( var.ne.'INPOLYGON' ) then
|
|
|
533 |
write(*,'(i5,f15.4)') i,out(i)
|
|
|
534 |
else
|
|
|
535 |
write(*,'(i5,f15.4,a2)') i,out(i),char(nint(out(i)))
|
|
|
536 |
endif
|
|
|
537 |
enddo
|
|
|
538 |
|
|
|
539 |
c Get times
|
|
|
540 |
i=i+1
|
|
|
541 |
ntim=nint(out(i))
|
|
|
542 |
write(*,'(i5,f15.4,7x,7x,a6)') i,out(i),'#TIMES'
|
|
|
543 |
do j=1,ntim
|
|
|
544 |
i=i+1
|
|
|
545 |
write(*,'(i5,f15.4,f7.0)') i,out(i),times(nint(out(i)))
|
|
|
546 |
enddo
|
|
|
547 |
|
|
|
548 |
c Get time mode
|
|
|
549 |
i=i+1
|
|
|
550 |
itime=nint(out(i))
|
|
|
551 |
if (itime.eq.1) strtim='ALL'
|
|
|
552 |
if (itime.eq.2) strtim='ANY'
|
|
|
553 |
if (itime.eq.3) strtim='NONE'
|
|
|
554 |
if (itime.lt.0) strtim='TRIGGER'
|
|
|
555 |
|
|
|
556 |
if ( strtim.ne.'TRIGGER' ) then
|
|
|
557 |
write(*,'(i5,f15.4,10x,a10)') i,out(i),trim(strtim)
|
|
|
558 |
else
|
|
|
559 |
write(*,'(i5,f15.4,10x,a10,a3,i3)') i,out(i),trim(strtim),
|
|
|
560 |
> ' ->',abs(itime)
|
|
|
561 |
endif
|
|
|
562 |
|
|
|
563 |
goto 100
|
|
|
564 |
|
|
|
565 |
endif
|
|
|
566 |
|
|
|
567 |
c Exit point
|
|
|
568 |
200 continue
|
|
|
569 |
|
|
|
570 |
write(*,*) '---------------------------------------'
|
|
|
571 |
|
|
|
572 |
end
|
|
|
573 |
|
|
|
574 |
|
|
|
575 |
c --------------------------------------------------------------
|
|
|
576 |
c Read and decode a selection set
|
|
|
577 |
c --------------------------------------------------------------
|
|
|
578 |
|
|
|
579 |
subroutine decode(fid,out,n,vars,nvars,times,ntimes,regionf)
|
|
|
580 |
|
|
|
581 |
c A selection file is opened with file id <fid> and transformed
|
|
|
582 |
c into a set of commands applied to the trajectories. On input
|
|
|
583 |
c <n> sets the maximum dimension of <out>, on output it gives the
|
|
|
584 |
c total length of the command string. The output is a list of
|
|
|
585 |
c commands with the following format:
|
|
|
586 |
c
|
|
|
587 |
c out(i) = Command
|
|
|
588 |
c out(i+1) = Column index of variable
|
|
|
589 |
c out(i+2) = Mode for variable
|
|
|
590 |
c out(i+3) = Number of parameters (n)
|
|
|
591 |
c out(i+4:i+4+n) = Parameters
|
|
|
592 |
c out(i+5+n) = Number of times
|
|
|
593 |
c out(i+6+n:i+6+n+m) = List of time indices (m)
|
|
|
594 |
c out(i+7+n+m) = Time mode
|
|
|
595 |
c
|
|
|
596 |
c For SPECIAL commands (to be coded in <special.f>) the format is:
|
|
|
597 |
c
|
|
|
598 |
c out(i) = Length of command string (n)
|
|
|
599 |
c out(i+1:i+n) = Command string
|
|
|
600 |
c out(i+n+1) = Number of parameters (m)
|
|
|
601 |
c out(i+n+2:i+n+1+m) = List of parameters
|
|
|
602 |
c
|
|
|
603 |
c The following coding is used
|
|
|
604 |
c
|
|
|
605 |
c Command Variable mode Time mode
|
|
|
606 |
c ---------- ------------- ---------
|
|
|
607 |
c GT 1 VALUE 1 ALL 1
|
|
|
608 |
c LT 2 MEAN 2 ANY 2
|
|
|
609 |
c IN 3 MAX 3 NONE 3
|
|
|
610 |
c OUT 4 MIN 4 TRIGGER -i (i the trigger index)
|
|
|
611 |
c EQ 5 VAR 5
|
|
|
612 |
c BEGIN -1 SUM 6
|
|
|
613 |
c END -2 CHANGE 7
|
|
|
614 |
c AND -3 DIFF 8
|
|
|
615 |
c OR -4
|
|
|
616 |
c TRUE 6
|
|
|
617 |
c FALSE 7
|
|
|
618 |
c SPECIAL 0
|
|
|
619 |
c ALL 8 (TRIGGER)
|
|
|
620 |
c ANY 9 (TRIGGER)
|
|
|
621 |
c NONE 10 (TRIGGER)
|
|
|
622 |
|
|
|
623 |
|
|
|
624 |
c Several "implicit variables" are supported - out(i+1):
|
|
|
625 |
c
|
|
|
626 |
c DIST -1 : Path length of the trajectory
|
|
|
627 |
c DIST0 -2 : Distance from starting position
|
|
|
628 |
c INPOLYGON -3 : Specified polygon region
|
|
|
629 |
c INBOX -4 : Longitude/latitude rectangle
|
|
|
630 |
c INCIRCLE -5 : Within a specified radius
|
|
|
631 |
c INREGION -6 : Within a specified rehion on the region file
|
|
|
632 |
c TRIGGER -7 : Trigger field
|
|
|
633 |
c
|
|
|
634 |
c For the special commands BEGIN, END, AND and OR, only one field
|
|
|
635 |
c in <out> is used.
|
|
|
636 |
|
|
|
637 |
implicit none
|
|
|
638 |
|
|
|
639 |
c Declaration of subroutine parameters
|
|
|
640 |
integer fid
|
|
|
641 |
integer n
|
|
|
642 |
real out(n)
|
|
|
643 |
integer nvars
|
|
|
644 |
character*80 vars(nvars)
|
|
|
645 |
integer ntimes
|
|
|
646 |
real times(ntimes)
|
|
|
647 |
character*80 regionf
|
|
|
648 |
|
|
|
649 |
c Numerical epsilon
|
|
|
650 |
real eps
|
|
|
651 |
parameter (eps=0.001)
|
|
|
652 |
|
|
|
653 |
c A single command term
|
|
|
654 |
character*80 cmd
|
|
|
655 |
character*80 var,mode
|
|
|
656 |
integer nval
|
|
|
657 |
real val(1000)
|
|
|
658 |
integer ntim
|
|
|
659 |
real tim(1000)
|
|
|
660 |
character*80 tmode
|
|
|
661 |
|
|
|
662 |
c Specification of a polygon
|
|
|
663 |
character*80 filename
|
|
|
664 |
integer pn ! Number of entries in lat/lon poly
|
|
|
665 |
real latpoly(1000) ! List of polygon latitudes
|
|
|
666 |
real lonpoly(1000) ! List of polygon longitudes
|
|
|
667 |
real loninpoly,latinpoly ! Lon/lat inside polygon
|
|
|
668 |
|
|
|
669 |
c Specification of a region
|
|
|
670 |
real xcorner(4)
|
|
|
671 |
real ycorner(4)
|
|
|
672 |
integer iregion
|
|
|
673 |
character*80 string
|
|
|
674 |
|
|
|
675 |
c Transformation to UPN (handling of logical operators)
|
|
|
676 |
integer nlogical
|
|
|
677 |
integer ilogical(n)
|
|
|
678 |
real tmp(n)
|
|
|
679 |
integer mlogical
|
|
|
680 |
integer isor
|
|
|
681 |
|
|
|
682 |
c Auxiliary variables
|
|
|
683 |
integer i,j
|
|
|
684 |
integer j1,j2
|
|
|
685 |
integer flag(ntimes)
|
|
|
686 |
integer count
|
|
|
687 |
integer ok
|
|
|
688 |
integer itrigger
|
|
|
689 |
|
|
|
690 |
c Common block for initialisation of polygon check
|
|
|
691 |
real tlonv(1000),vlat_c(1000),vlon_c(1000)
|
|
|
692 |
real xlat_c,xlon_c
|
|
|
693 |
integer ibndry,nv_c
|
|
|
694 |
common /spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry
|
|
|
695 |
|
|
|
696 |
c ------ Decode single commands ---------------------
|
|
|
697 |
|
|
|
698 |
c Reset the filename for polygons
|
|
|
699 |
filename='nil'
|
|
|
700 |
|
|
|
701 |
c Reset the counter for logical commands
|
|
|
702 |
nlogical=0
|
|
|
703 |
|
|
|
704 |
c Loop through all commands
|
|
|
705 |
100 continue
|
|
|
706 |
|
|
|
707 |
c Read next command (handle special cases)
|
|
|
708 |
read(fid,*) cmd
|
|
|
709 |
|
|
|
710 |
c Special handling of SPECIAL commands
|
|
|
711 |
if ( cmd.eq.'SPECIAL' ) then
|
|
|
712 |
|
|
|
713 |
c Set the flag for SPECIAL command
|
|
|
714 |
n = 1
|
|
|
715 |
out(n) = 0.
|
|
|
716 |
|
|
|
717 |
c Read the command string
|
|
|
718 |
read(fid,*) var
|
|
|
719 |
|
|
|
720 |
c Add the command string
|
|
|
721 |
n = n + 1
|
|
|
722 |
out(n) = len_trim(var)
|
|
|
723 |
do j=1,len_trim(var)
|
|
|
724 |
n = n + 1
|
|
|
725 |
out(n) = ichar(var(j:j))
|
|
|
726 |
enddo
|
|
|
727 |
|
|
|
728 |
c Read the parameters
|
|
|
729 |
read(fid,*) nval
|
|
|
730 |
read(fid,*) (val(i),i=1,nval)
|
|
|
731 |
|
|
|
732 |
c Add the parameters
|
|
|
733 |
n = n + 1
|
|
|
734 |
out(n)=real(nval)
|
|
|
735 |
do i=1,nval
|
|
|
736 |
n=n+1
|
|
|
737 |
out(n)=val(i)
|
|
|
738 |
enddo
|
|
|
739 |
|
|
|
740 |
c Goto exit point - nothing more top do
|
|
|
741 |
goto 350
|
|
|
742 |
|
|
|
743 |
endif
|
|
|
744 |
|
|
|
745 |
c Handle structure commands
|
|
|
746 |
if ( cmd.eq.'BEGIN') then
|
|
|
747 |
out(1)=-1.
|
|
|
748 |
n=1
|
|
|
749 |
nlogical=1
|
|
|
750 |
ilogical(1)=n
|
|
|
751 |
goto 200
|
|
|
752 |
elseif ( cmd.eq.'AND' ) then
|
|
|
753 |
n=n+1
|
|
|
754 |
out(n)=-3.
|
|
|
755 |
nlogical=nlogical+1
|
|
|
756 |
ilogical(nlogical)=n
|
|
|
757 |
goto 200
|
|
|
758 |
elseif ( cmd.eq.'OR' ) then
|
|
|
759 |
n=n+1
|
|
|
760 |
out(n)=-4.
|
|
|
761 |
nlogical=nlogical+1
|
|
|
762 |
ilogical(nlogical)=n
|
|
|
763 |
goto 200
|
|
|
764 |
elseif ( cmd.eq.'END' ) then
|
|
|
765 |
n=n+1
|
|
|
766 |
out(n)=-2.
|
|
|
767 |
nlogical=nlogical+1
|
|
|
768 |
ilogical(nlogical)=n
|
|
|
769 |
goto 300
|
|
|
770 |
endif
|
|
|
771 |
|
|
|
772 |
c Read other fields associated with the command
|
|
|
773 |
read(fid,*) var,mode
|
|
|
774 |
|
|
|
775 |
c Read parameter
|
|
|
776 |
if ( var.eq.'INPOLYGON' ) then
|
|
|
777 |
read(fid,*) nval
|
|
|
778 |
read(fid,*) filename
|
|
|
779 |
filename = trim(filename)
|
|
|
780 |
else
|
|
|
781 |
read(fid,*) nval
|
|
|
782 |
read(fid,*) (val(i),i=1,nval)
|
|
|
783 |
endif
|
|
|
784 |
|
|
|
785 |
c Read times
|
|
|
786 |
read(fid,*) ntim
|
|
|
787 |
read(fid,*) (tim(i),i=1,ntim)
|
|
|
788 |
read(fid,*) tmode
|
|
|
789 |
|
|
|
790 |
c Bring CAPITAL "TIME,LAT,LON,P" into "time,lat,lon,p"
|
|
|
791 |
if (var.eq.'TIME') var='time'
|
|
|
792 |
if (var.eq.'LAT' ) var='lat'
|
|
|
793 |
if (var.eq.'LON' ) var='lon'
|
|
|
794 |
if (var.eq.'P' ) var='p'
|
|
|
795 |
|
|
|
796 |
c If the time mode is 'TRIGGER', all times of a trajectory
|
|
|
797 |
c must be considered
|
|
|
798 |
if ( tmode.eq.'TRIGGER' ) then
|
|
|
799 |
itrigger = nint(tim(1))
|
|
|
800 |
ntim = 1
|
|
|
801 |
tim(1) = -999.
|
|
|
802 |
endif
|
|
|
803 |
|
|
|
804 |
c Special times: transform into real time
|
|
|
805 |
do i=1,ntim
|
|
|
806 |
if ( abs(tim(i)+996.).lt.eps ) tim(i)=times(1)
|
|
|
807 |
if ( abs(tim(i)+995.).lt.eps ) tim(i)=times(ntimes)
|
|
|
808 |
enddo
|
|
|
809 |
|
|
|
810 |
c Check whether times are valid
|
|
|
811 |
ok=0
|
|
|
812 |
do i=1,ntim
|
|
|
813 |
if ( (abs(tim(i)+994.).gt.eps).and.
|
|
|
814 |
> (abs(tim(i)+999.).gt.eps) )
|
|
|
815 |
> then
|
|
|
816 |
do j=1,ntimes
|
|
|
817 |
if ( abs(tim(i)-times(j)).lt.eps ) then
|
|
|
818 |
ok=ok+1
|
|
|
819 |
endif
|
|
|
820 |
enddo
|
|
|
821 |
else
|
|
|
822 |
ok=ok+1
|
|
|
823 |
endif
|
|
|
824 |
enddo
|
|
|
825 |
if (ok.ne.ntim) goto 400
|
|
|
826 |
|
|
|
827 |
c Select all times which are included in the criterion
|
|
|
828 |
do i=1,ntimes
|
|
|
829 |
flag(i)=0
|
|
|
830 |
enddo
|
|
|
831 |
i=1
|
|
|
832 |
150 if (i.le.ntim) then
|
|
|
833 |
|
|
|
834 |
c A list of times
|
|
|
835 |
if ( (abs(tim(i)+994.).lt.eps) ) then
|
|
|
836 |
j1=0
|
|
|
837 |
do j=1,ntimes
|
|
|
838 |
if ( abs(tim(i-1)-times(j)).lt.eps ) then
|
|
|
839 |
j1=j
|
|
|
840 |
endif
|
|
|
841 |
enddo
|
|
|
842 |
j2=0
|
|
|
843 |
do j=1,ntimes
|
|
|
844 |
if ( abs(tim(i+1)-times(j)).lt.eps ) then
|
|
|
845 |
j2=j
|
|
|
846 |
endif
|
|
|
847 |
enddo
|
|
|
848 |
if ( (j1.eq.0).or.(j2.eq.0) ) goto 400
|
|
|
849 |
do j=j1,j2
|
|
|
850 |
flag(j)=1
|
|
|
851 |
enddo
|
|
|
852 |
i=i+1
|
|
|
853 |
|
|
|
854 |
c Explicitly given time value
|
|
|
855 |
else
|
|
|
856 |
do j=1,ntimes
|
|
|
857 |
if ( abs(tim(i)-times(j)).lt.eps ) then
|
|
|
858 |
flag(j)=1
|
|
|
859 |
endif
|
|
|
860 |
enddo
|
|
|
861 |
|
|
|
862 |
endif
|
|
|
863 |
|
|
|
864 |
i=i+1
|
|
|
865 |
goto 150
|
|
|
866 |
|
|
|
867 |
endif
|
|
|
868 |
|
|
|
869 |
c Write command identifier
|
|
|
870 |
n=n+1
|
|
|
871 |
if (cmd.eq.'GT' ) out(n)= 1.
|
|
|
872 |
if (cmd.eq.'LT' ) out(n)= 2.
|
|
|
873 |
if (cmd.eq.'IN' ) out(n)= 3.
|
|
|
874 |
if (cmd.eq.'OUT' ) out(n)= 4.
|
|
|
875 |
if (cmd.eq.'EQ' ) out(n)= 5.
|
|
|
876 |
if (cmd.eq.'TRUE' ) out(n)= 6.
|
|
|
877 |
if (cmd.eq.'FALSE' ) out(n)= 7.
|
|
|
878 |
if (cmd.eq.'ALL ' ) out(n)= 8.
|
|
|
879 |
if (cmd.eq.'ANY ' ) out(n)= 9.
|
|
|
880 |
if (cmd.eq.'NONE ' ) out(n)=10.
|
|
|
881 |
|
|
|
882 |
c Write index for variable - force implicit trigger
|
|
|
883 |
ok=0
|
|
|
884 |
do j=1,nvars
|
|
|
885 |
if (vars(j).eq.var) ok=j
|
|
|
886 |
enddo
|
|
|
887 |
|
|
|
888 |
if (var.eq.'TRIGGER') ok = 0
|
|
|
889 |
|
|
|
890 |
if (ok.eq.0) then
|
|
|
891 |
if (var.eq.'DIST') then
|
|
|
892 |
ok = -1
|
|
|
893 |
elseif (var.eq.'DIST0') then
|
|
|
894 |
ok = -2
|
|
|
895 |
elseif (var.eq.'INPOLYGON') then
|
|
|
896 |
ok = -3
|
|
|
897 |
elseif (var.eq.'INBOX') then
|
|
|
898 |
ok = -4
|
|
|
899 |
elseif (var.eq.'INCIRCLE') then
|
|
|
900 |
ok = -5
|
|
|
901 |
elseif (var.eq.'INREGION') then
|
|
|
902 |
ok = -6
|
|
|
903 |
elseif (var.eq.'TRIGGER') then
|
|
|
904 |
ok = -7
|
|
|
905 |
else
|
|
|
906 |
goto 400
|
|
|
907 |
endif
|
|
|
908 |
endif
|
|
|
909 |
n=n+1
|
|
|
910 |
out(n)=real(ok)
|
|
|
911 |
|
|
|
912 |
c Write mode for variable
|
|
|
913 |
ok=0
|
|
|
914 |
if (mode.eq.'VALUE' ) ok=1
|
|
|
915 |
if (mode.eq.'MEAN' ) ok=2
|
|
|
916 |
if (mode.eq.'MAX' ) ok=3
|
|
|
917 |
if (mode.eq.'MIN' ) ok=4
|
|
|
918 |
if (mode.eq.'VAR' ) ok=5
|
|
|
919 |
if (mode.eq.'SUM' ) ok=6
|
|
|
920 |
if (mode.eq.'CHANGE' ) ok=7
|
|
|
921 |
if (mode.eq.'DIFF' ) ok=8
|
|
|
922 |
if (ok.eq.0) goto 400
|
|
|
923 |
n=n+1
|
|
|
924 |
out(n)=real(ok)
|
|
|
925 |
|
|
|
926 |
c Write the parameter values: INPOLYGON
|
|
|
927 |
if ( var.eq.'INPOLYGON' ) then
|
|
|
928 |
|
|
|
929 |
n = n+1
|
|
|
930 |
out(n) = len_trim(filename)
|
|
|
931 |
do j=1,len_trim(filename)
|
|
|
932 |
n = n + 1
|
|
|
933 |
out(n) = ichar(filename(j:j))
|
|
|
934 |
enddo
|
|
|
935 |
|
|
|
936 |
c Write parameter value: INREGION
|
|
|
937 |
elseif ( var.eq.'INREGION' ) then
|
|
|
938 |
|
|
|
939 |
iregion = nint(val(1))
|
|
|
940 |
|
|
|
941 |
open(fid+1,file=regionf)
|
|
|
942 |
|
|
|
943 |
50 read(fid+1,*,end=51) string
|
|
|
944 |
|
|
|
945 |
if ( string(1:1).ne.'#' ) then
|
|
|
946 |
call regionsplit(string,i,xcorner,ycorner)
|
|
|
947 |
if ( i.eq.iregion ) goto 52
|
|
|
948 |
endif
|
|
|
949 |
|
|
|
950 |
goto 50
|
|
|
951 |
|
|
|
952 |
51 close(fid+1)
|
|
|
953 |
|
|
|
954 |
print*,' ERROR: region ',iregion,' not found on ',
|
|
|
955 |
> trim(regionf)
|
|
|
956 |
stop
|
|
|
957 |
|
|
|
958 |
52 continue
|
|
|
959 |
|
|
|
960 |
n = n + 1
|
|
|
961 |
out(n) = 8 ! Number of parameters
|
|
|
962 |
do i=1,4
|
|
|
963 |
n = n + 1
|
|
|
964 |
out(n) = xcorner(i)
|
|
|
965 |
enddo
|
|
|
966 |
do i=1,4
|
|
|
967 |
n = n + 1
|
|
|
968 |
out(n) = ycorner(i)
|
|
|
969 |
enddo
|
|
|
970 |
|
|
|
971 |
c Write parameter values: all other cases
|
|
|
972 |
else
|
|
|
973 |
n=n+1
|
|
|
974 |
out(n)=real(nval)
|
|
|
975 |
do i=1,nval
|
|
|
976 |
n=n+1
|
|
|
977 |
out(n)=val(i)
|
|
|
978 |
enddo
|
|
|
979 |
endif
|
|
|
980 |
|
|
|
981 |
c All times are selected
|
|
|
982 |
if ( abs(tim(1)+999.).lt.eps ) then
|
|
|
983 |
n=n+1
|
|
|
984 |
out(n)=real(ntimes)
|
|
|
985 |
do i=1,ntimes
|
|
|
986 |
n=n+1
|
|
|
987 |
out(n)=real(i)
|
|
|
988 |
enddo
|
|
|
989 |
|
|
|
990 |
c A selection of times is given
|
|
|
991 |
else
|
|
|
992 |
count=0
|
|
|
993 |
do i=1,ntimes
|
|
|
994 |
count=count+flag(i)
|
|
|
995 |
enddo
|
|
|
996 |
n=n+1
|
|
|
997 |
out(n)=real(count)
|
|
|
998 |
do i=1,ntimes
|
|
|
999 |
if (flag(i).eq.1) then
|
|
|
1000 |
n=n+1
|
|
|
1001 |
out(n)=real(i)
|
|
|
1002 |
endif
|
|
|
1003 |
enddo
|
|
|
1004 |
endif
|
|
|
1005 |
|
|
|
1006 |
c Write the time mode
|
|
|
1007 |
if ( tmode.eq.'ALL') then
|
|
|
1008 |
n=n+1
|
|
|
1009 |
out(n)=1.
|
|
|
1010 |
elseif ( tmode.eq.'ANY') then
|
|
|
1011 |
n=n+1
|
|
|
1012 |
out(n)=2.
|
|
|
1013 |
elseif ( tmode.eq.'NONE') then
|
|
|
1014 |
n=n+1
|
|
|
1015 |
out(n)=3.
|
|
|
1016 |
elseif ( tmode.eq.'TRIGGER') then
|
|
|
1017 |
n=n+1
|
|
|
1018 |
out(n)=-real(itrigger)
|
|
|
1019 |
endif
|
|
|
1020 |
|
|
|
1021 |
c End loop: handle single command
|
|
|
1022 |
200 continue
|
|
|
1023 |
goto 100
|
|
|
1024 |
|
|
|
1025 |
c End loop: loop over all commands
|
|
|
1026 |
300 continue
|
|
|
1027 |
|
|
|
1028 |
c ------ Read polygon file, if requested -----------
|
|
|
1029 |
if ( filename.ne.'nil' ) then
|
|
|
1030 |
|
|
|
1031 |
print*
|
|
|
1032 |
print*,
|
|
|
1033 |
> '---- POLYGON --------------------------------------------'
|
|
|
1034 |
|
|
|
1035 |
print*
|
|
|
1036 |
print*,'Filename = ',trim(filename)
|
|
|
1037 |
print*
|
|
|
1038 |
|
|
|
1039 |
c Read list of polygon coordinates from file
|
|
|
1040 |
pn = 0
|
|
|
1041 |
open(fid+1,file=filename)
|
|
|
1042 |
read(fid+1,*) loninpoly,latinpoly
|
|
|
1043 |
print*,'Inside (lon/lat) =',loninpoly,latinpoly
|
|
|
1044 |
print*
|
|
|
1045 |
510 continue
|
|
|
1046 |
pn = pn + 1
|
|
|
1047 |
read(fid+1,*,end=511) lonpoly(pn),
|
|
|
1048 |
> latpoly(pn)
|
|
|
1049 |
|
|
|
1050 |
print*,pn,lonpoly(pn),latpoly(pn)
|
|
|
1051 |
|
|
|
1052 |
goto 510
|
|
|
1053 |
511 continue
|
|
|
1054 |
pn = pn - 1
|
|
|
1055 |
close(fid+1)
|
|
|
1056 |
|
|
|
1057 |
c Define the polygon boundaries
|
|
|
1058 |
call DefSPolyBndry(latpoly,lonpoly,pn,latinpoly,loninpoly)
|
|
|
1059 |
|
|
|
1060 |
endif
|
|
|
1061 |
|
|
|
1062 |
c ------ Transform to UPN --------------------------
|
|
|
1063 |
|
|
|
1064 |
c Check whether logical commands are ok
|
|
|
1065 |
mlogical=nint(out(ilogical(1)))
|
|
|
1066 |
if ( mlogical.ne.-1) goto 400
|
|
|
1067 |
mlogical=nint(out(ilogical(nlogical)))
|
|
|
1068 |
if ( mlogical.ne.-2) goto 400
|
|
|
1069 |
|
|
|
1070 |
c No transformation necessary if only one command
|
|
|
1071 |
if (nlogical.eq.2) goto 350
|
|
|
1072 |
|
|
|
1073 |
c Copy the output to temporary list
|
|
|
1074 |
do i=1,n
|
|
|
1075 |
tmp(i)=out(i)
|
|
|
1076 |
enddo
|
|
|
1077 |
|
|
|
1078 |
c Set BEGIN statement
|
|
|
1079 |
n=1
|
|
|
1080 |
out(n)=-1.
|
|
|
1081 |
|
|
|
1082 |
c Reorder commands and transform to UPN
|
|
|
1083 |
isor=0
|
|
|
1084 |
do i=1,nlogical-1
|
|
|
1085 |
|
|
|
1086 |
c Get the logical command
|
|
|
1087 |
mlogical=nint(out(ilogical(i)))
|
|
|
1088 |
|
|
|
1089 |
c Connecting OR
|
|
|
1090 |
if (mlogical.eq.-4) then
|
|
|
1091 |
if (isor.eq.1) then
|
|
|
1092 |
n=n+1
|
|
|
1093 |
out(n)=-4.
|
|
|
1094 |
else
|
|
|
1095 |
isor=1
|
|
|
1096 |
endif
|
|
|
1097 |
endif
|
|
|
1098 |
|
|
|
1099 |
c Put the command onto the stack
|
|
|
1100 |
do j=ilogical(i)+1,ilogical(i+1)-1
|
|
|
1101 |
n=n+1
|
|
|
1102 |
out(n)=tmp(j)
|
|
|
1103 |
enddo
|
|
|
1104 |
|
|
|
1105 |
c Connecting AND operator
|
|
|
1106 |
if ( mlogical.eq.-3 ) then
|
|
|
1107 |
n=n+1
|
|
|
1108 |
out(n)=-3.
|
|
|
1109 |
endif
|
|
|
1110 |
|
|
|
1111 |
enddo
|
|
|
1112 |
|
|
|
1113 |
c Set final connecting OR
|
|
|
1114 |
if (isor.eq.1) then
|
|
|
1115 |
n=n+1
|
|
|
1116 |
out(n)=-4.
|
|
|
1117 |
endif
|
|
|
1118 |
|
|
|
1119 |
c Set END statement
|
|
|
1120 |
n=n+1
|
|
|
1121 |
out(n)=-2.
|
|
|
1122 |
|
|
|
1123 |
c ------ Exit point ---------------------------------
|
|
|
1124 |
|
|
|
1125 |
350 continue
|
|
|
1126 |
return
|
|
|
1127 |
|
|
|
1128 |
c ----- Exception handling --------------------------
|
|
|
1129 |
400 print*,'Invalid selection criterion... Stop'
|
|
|
1130 |
stop
|
|
|
1131 |
|
|
|
1132 |
end
|
|
|
1133 |
|
|
|
1134 |
|
|
|
1135 |
c --------------------------------------------------------------
|
|
|
1136 |
c Decide whether a trajectory is selected or not
|
|
|
1137 |
c --------------------------------------------------------------
|
|
|
1138 |
|
|
|
1139 |
subroutine select_tra (select,cmd,n,tra,trigger,ntim,ncol)
|
|
|
1140 |
|
|
|
1141 |
c Decide whether a single trajectory is selected (<select=1>) or
|
|
|
1142 |
c is not selected <select=0> according to the selection criterion
|
|
|
1143 |
c given in <cmd(ncmd)>. The selection criterion <cmd(ncmd)> is
|
|
|
1144 |
c returned from the call to the subroutine <decode>. The trajectory
|
|
|
1145 |
c is given in <tra(ntim,ncol)> where <ntim> is the number of times
|
|
|
1146 |
c and <ncol> is the number of columns.
|
|
|
1147 |
c
|
|
|
1148 |
c Important note: the structure of <tra(ntim,ncol)> must match to the
|
|
|
1149 |
c call parameter <vars,nvars,times,ntimes> in subroutine <decode>.
|
|
|
1150 |
|
|
|
1151 |
implicit none
|
|
|
1152 |
|
|
|
1153 |
c Declaration of subroutine parameters
|
|
|
1154 |
integer select
|
|
|
1155 |
integer n
|
|
|
1156 |
real cmd(n)
|
|
|
1157 |
integer ntim,ncol
|
|
|
1158 |
real tra(ntim,ncol)
|
|
|
1159 |
integer trigger(ntim)
|
|
|
1160 |
|
|
|
1161 |
c Numerical epsilon (for test of equality)
|
|
|
1162 |
real eps
|
|
|
1163 |
parameter (eps=0.000001)
|
|
|
1164 |
|
|
|
1165 |
c A single command and the associated field
|
|
|
1166 |
integer icmd,ivar,imode,itime,nsel,nval
|
|
|
1167 |
integer time(ntim)
|
|
|
1168 |
real param(100)
|
|
|
1169 |
real var (ntim)
|
|
|
1170 |
integer intvar(ntim)
|
|
|
1171 |
|
|
|
1172 |
c Boolean values for a single time, a single command and build-up
|
|
|
1173 |
integer stack(100)
|
|
|
1174 |
integer nstack
|
|
|
1175 |
integer istrue(ntim)
|
|
|
1176 |
integer decision
|
|
|
1177 |
|
|
|
1178 |
c Auxiliary variables
|
|
|
1179 |
integer i,j,k
|
|
|
1180 |
real tmp,mea
|
|
|
1181 |
integer istack1,istack2
|
|
|
1182 |
real lat0,lon0,lat1,lon1
|
|
|
1183 |
real length(ntim)
|
|
|
1184 |
integer flag
|
|
|
1185 |
real dist
|
|
|
1186 |
real xcorner(4),ycorner(4)
|
|
|
1187 |
integer iparam
|
|
|
1188 |
character ch
|
|
|
1189 |
|
|
|
1190 |
c Common block for initialisation of polygon check
|
|
|
1191 |
real tlonv(1000),vlat_c(1000),vlon_c(1000)
|
|
|
1192 |
real xlat_c,xlon_c
|
|
|
1193 |
integer ibndry,nv_c
|
|
|
1194 |
common /spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry
|
|
|
1195 |
|
|
|
1196 |
c Externals
|
|
|
1197 |
real sdis ! Spherical distance
|
|
|
1198 |
external sdis
|
|
|
1199 |
integer inregion ! Boolean flag for regions
|
|
|
1200 |
external inregion
|
|
|
1201 |
|
|
|
1202 |
c Reset the decision stack (with locical values)
|
|
|
1203 |
nstack=0
|
|
|
1204 |
|
|
|
1205 |
c Loop through the complete command list
|
|
|
1206 |
i=0
|
|
|
1207 |
100 if (i.lt.n) then
|
|
|
1208 |
|
|
|
1209 |
c --- Get the command -------------------------------
|
|
|
1210 |
i=i+1
|
|
|
1211 |
icmd=nint(cmd(i))
|
|
|
1212 |
|
|
|
1213 |
c --- Handle structural commands (BEGIN, END, AND, OR)
|
|
|
1214 |
|
|
|
1215 |
c Handle BEGIN statement
|
|
|
1216 |
if ( icmd.eq.-1) then
|
|
|
1217 |
nstack=0
|
|
|
1218 |
goto 200
|
|
|
1219 |
endif
|
|
|
1220 |
|
|
|
1221 |
c Handle END statement
|
|
|
1222 |
if (icmd.eq.-2) then
|
|
|
1223 |
goto 300
|
|
|
1224 |
endif
|
|
|
1225 |
|
|
|
1226 |
c Handle AND statement
|
|
|
1227 |
if (icmd.eq.-3) then
|
|
|
1228 |
istack1=stack(nstack)
|
|
|
1229 |
nstack=nstack-1
|
|
|
1230 |
istack2=stack(nstack)
|
|
|
1231 |
if ((istack1.eq.1).and.(istack2.eq.1)) then
|
|
|
1232 |
stack(nstack)=1
|
|
|
1233 |
else
|
|
|
1234 |
stack(nstack)=0
|
|
|
1235 |
endif
|
|
|
1236 |
goto 200
|
|
|
1237 |
endif
|
|
|
1238 |
|
|
|
1239 |
c Handle OR statement
|
|
|
1240 |
if (icmd.eq.-4) then
|
|
|
1241 |
istack1=stack(nstack)
|
|
|
1242 |
nstack=nstack-1
|
|
|
1243 |
istack2=stack(nstack)
|
|
|
1244 |
if ((istack1.eq.1).or.(istack2.eq.1)) then
|
|
|
1245 |
stack(nstack)=1
|
|
|
1246 |
else
|
|
|
1247 |
stack(nstack)=0
|
|
|
1248 |
endif
|
|
|
1249 |
goto 200
|
|
|
1250 |
endif
|
|
|
1251 |
|
|
|
1252 |
c --- Get all command details (parameters, modes, times)
|
|
|
1253 |
|
|
|
1254 |
c Get variable (<ivar> gets the column index in <tra>)
|
|
|
1255 |
i=i+1
|
|
|
1256 |
ivar=nint(cmd(i))
|
|
|
1257 |
|
|
|
1258 |
c Get variable mode
|
|
|
1259 |
i=i+1
|
|
|
1260 |
imode=nint(cmd(i))
|
|
|
1261 |
|
|
|
1262 |
c Get parameter values
|
|
|
1263 |
i=i+1
|
|
|
1264 |
nval=nint(cmd(i))
|
|
|
1265 |
do j=1,nval
|
|
|
1266 |
i=i+1
|
|
|
1267 |
param(j)=cmd(i)
|
|
|
1268 |
enddo
|
|
|
1269 |
|
|
|
1270 |
c Get times (<time(j)> gets the row indices of <tra>)
|
|
|
1271 |
i=i+1
|
|
|
1272 |
nsel=nint(cmd(i))
|
|
|
1273 |
do j=1,nsel
|
|
|
1274 |
i=i+1
|
|
|
1275 |
time(j)=nint(cmd(i))
|
|
|
1276 |
enddo
|
|
|
1277 |
|
|
|
1278 |
c Get time mode
|
|
|
1279 |
i=i+1
|
|
|
1280 |
itime=nint(cmd(i))
|
|
|
1281 |
|
|
|
1282 |
c --- Prepare field values for analysis -----------
|
|
|
1283 |
|
|
|
1284 |
c Implicit variable: DIST
|
|
|
1285 |
if ( ivar.eq. -1 ) then
|
|
|
1286 |
length(1) = 0.
|
|
|
1287 |
do j=2,ntim
|
|
|
1288 |
lon0 = tra(j-1,2)
|
|
|
1289 |
lat0 = tra(j-1,3)
|
|
|
1290 |
lon1 = tra(j ,2)
|
|
|
1291 |
lat1 = tra(j ,3)
|
|
|
1292 |
length(j) = length(j-1) + sdis(lon0,lat0,lon1,lat1)
|
|
|
1293 |
enddo
|
|
|
1294 |
do j=1,nsel
|
|
|
1295 |
var(j) = length( time(j) )
|
|
|
1296 |
enddo
|
|
|
1297 |
|
|
|
1298 |
|
|
|
1299 |
c Implict variable: DIST0
|
|
|
1300 |
elseif ( ivar.eq. -2 ) then
|
|
|
1301 |
do j=1,nsel
|
|
|
1302 |
lon0 = tra(1 ,2)
|
|
|
1303 |
lat0 = tra(1 ,3)
|
|
|
1304 |
lon1 = tra(time(j),2)
|
|
|
1305 |
lat1 = tra(time(j),3)
|
|
|
1306 |
var(j) = sdis(lon0,lat0,lon1,lat1)
|
|
|
1307 |
enddo
|
|
|
1308 |
|
|
|
1309 |
c Implict variable: INPOLYGON
|
|
|
1310 |
elseif ( ivar.eq. -3 ) then
|
|
|
1311 |
do j=1,nsel
|
|
|
1312 |
lon1 = tra(time(j),2)
|
|
|
1313 |
lat1 = tra(time(j),3)
|
|
|
1314 |
call LctPtRelBndry(lat1,lon1,flag)
|
|
|
1315 |
if ( (flag.eq.1).or.(flag.eq.2) ) then
|
|
|
1316 |
var(j) = 1.
|
|
|
1317 |
else
|
|
|
1318 |
var(j) = 0.
|
|
|
1319 |
endif
|
|
|
1320 |
enddo
|
|
|
1321 |
|
|
|
1322 |
c Implict variable: INBOX
|
|
|
1323 |
elseif ( ivar.eq. -4 ) then
|
|
|
1324 |
do j=1,nsel
|
|
|
1325 |
lon1 = tra(time(j),2)
|
|
|
1326 |
lat1 = tra(time(j),3)
|
|
|
1327 |
if ( ( lon1.ge.param(1) ).and. ! lonmin
|
|
|
1328 |
> ( lon1.le.param(2) ).and. ! lonmax
|
|
|
1329 |
> ( lat1.ge.param(3) ).and. ! latmin
|
|
|
1330 |
> ( lat1.le.param(4) ) ) ! latmax
|
|
|
1331 |
> then
|
|
|
1332 |
var(j) = 1
|
|
|
1333 |
else
|
|
|
1334 |
var(j) = 0
|
|
|
1335 |
endif
|
|
|
1336 |
enddo
|
|
|
1337 |
|
|
|
1338 |
c Implict variable: INCIRCLE (lonc=param(1),latc=param(2),radius=param(3))
|
|
|
1339 |
elseif ( ivar.eq. -5 ) then
|
|
|
1340 |
do j=1,nsel
|
|
|
1341 |
lon1 = tra(time(j),2)
|
|
|
1342 |
lat1 = tra(time(j),3)
|
|
|
1343 |
|
|
|
1344 |
dist = sdis( lon1,lat1,param(1),param(2) )
|
|
|
1345 |
|
|
|
1346 |
if ( dist.le.param(3) ) then
|
|
|
1347 |
var(j) = 1
|
|
|
1348 |
else
|
|
|
1349 |
var(j) = 0
|
|
|
1350 |
endif
|
|
|
1351 |
enddo
|
|
|
1352 |
|
|
|
1353 |
c Implict variable: INREGION (xcorner=param(1..4),ycorner=param(5..8) )
|
|
|
1354 |
elseif ( ivar.eq.-6 ) then
|
|
|
1355 |
|
|
|
1356 |
do j=1,4
|
|
|
1357 |
xcorner(j) = param(j )
|
|
|
1358 |
ycorner(j) = param(j+4)
|
|
|
1359 |
enddo
|
|
|
1360 |
|
|
|
1361 |
do j=1,nsel
|
|
|
1362 |
lon1 = tra(time(j),2)
|
|
|
1363 |
lat1 = tra(time(j),3)
|
|
|
1364 |
var(j) = inregion (lon1,lat1,xcorner,ycorner)
|
|
|
1365 |
enddo
|
|
|
1366 |
|
|
|
1367 |
c Implict variable: TRIGGER
|
|
|
1368 |
elseif ( ivar.eq. -7 ) then
|
|
|
1369 |
do j=1,nsel
|
|
|
1370 |
intvar(j) = trigger( time(j) )
|
|
|
1371 |
enddo
|
|
|
1372 |
|
|
|
1373 |
|
|
|
1374 |
c Explicit variable (column index <ivar>)
|
|
|
1375 |
else
|
|
|
1376 |
do j=1,nsel
|
|
|
1377 |
var(j) = tra(time(j),ivar)
|
|
|
1378 |
enddo
|
|
|
1379 |
|
|
|
1380 |
endif
|
|
|
1381 |
|
|
|
1382 |
c Take MEAN of the variable (mean of selected times)
|
|
|
1383 |
if (imode.eq.2) then
|
|
|
1384 |
tmp=0.
|
|
|
1385 |
do j=1,nsel
|
|
|
1386 |
tmp=tmp+var(j)
|
|
|
1387 |
enddo
|
|
|
1388 |
var(1)=tmp/real(nsel)
|
|
|
1389 |
nsel=1
|
|
|
1390 |
|
|
|
1391 |
c Take MAX of the variable (maximum of selected times)
|
|
|
1392 |
elseif (imode.eq.3) then
|
|
|
1393 |
tmp=var(1)
|
|
|
1394 |
do j=2,nsel
|
|
|
1395 |
if (var(j).gt.tmp) tmp=var(j)
|
|
|
1396 |
enddo
|
|
|
1397 |
var(1)=tmp
|
|
|
1398 |
nsel=1
|
|
|
1399 |
|
|
|
1400 |
c Take MIN of the variable (minimum of selected times)
|
|
|
1401 |
elseif (imode.eq.4) then
|
|
|
1402 |
tmp=var(1)
|
|
|
1403 |
do j=2,nsel
|
|
|
1404 |
if (var(j).lt.tmp) tmp=var(j)
|
|
|
1405 |
enddo
|
|
|
1406 |
var(1)=tmp
|
|
|
1407 |
nsel=1
|
|
|
1408 |
|
|
|
1409 |
c Take VAR of the variable (variance over all selected times)
|
|
|
1410 |
elseif (imode.eq.5) then
|
|
|
1411 |
tmp=0.
|
|
|
1412 |
do j=1,nsel
|
|
|
1413 |
tmp=tmp+var(j)
|
|
|
1414 |
enddo
|
|
|
1415 |
mea=tmp/real(nsel)
|
|
|
1416 |
do j=1,nsel
|
|
|
1417 |
tmp=tmp+(var(j)-mea)**2
|
|
|
1418 |
enddo
|
|
|
1419 |
var(1)=1./real(nsel-1)*tmp
|
|
|
1420 |
nsel=1
|
|
|
1421 |
|
|
|
1422 |
c Take SUM of the variable (sum over all selected times)
|
|
|
1423 |
elseif (imode.eq.6) then
|
|
|
1424 |
tmp=0.
|
|
|
1425 |
do j=1,nsel
|
|
|
1426 |
tmp=tmp+var(j)
|
|
|
1427 |
enddo
|
|
|
1428 |
var(1)=tmp
|
|
|
1429 |
nsel=1
|
|
|
1430 |
|
|
|
1431 |
c Take CHANGE of the variable (change between first and last time)
|
|
|
1432 |
elseif (imode.eq.7) then
|
|
|
1433 |
var(1)=abs(var(1)-var(nsel))
|
|
|
1434 |
nsel=1
|
|
|
1435 |
|
|
|
1436 |
c Take DIFF of the variable (first minus last time)
|
|
|
1437 |
elseif (imode.eq.8) then
|
|
|
1438 |
var(1)=var(1)-var(nsel)
|
|
|
1439 |
nsel=1
|
|
|
1440 |
|
|
|
1441 |
endif
|
|
|
1442 |
|
|
|
1443 |
c --- Apply the operators to the single values ---
|
|
|
1444 |
|
|
|
1445 |
do j=1,nsel
|
|
|
1446 |
|
|
|
1447 |
c GT
|
|
|
1448 |
if (icmd.eq.1) then
|
|
|
1449 |
if (var(j).gt.param(1)) then
|
|
|
1450 |
istrue(j)=1
|
|
|
1451 |
else
|
|
|
1452 |
istrue(j)=0
|
|
|
1453 |
endif
|
|
|
1454 |
|
|
|
1455 |
c LT
|
|
|
1456 |
elseif (icmd.eq.2) then
|
|
|
1457 |
if (var(j).lt.param(1)) then
|
|
|
1458 |
istrue(j)=1
|
|
|
1459 |
else
|
|
|
1460 |
istrue(j)=0
|
|
|
1461 |
endif
|
|
|
1462 |
|
|
|
1463 |
c IN
|
|
|
1464 |
elseif (icmd.eq.3) then
|
|
|
1465 |
if ( (var(j).gt.param(1)).and.
|
|
|
1466 |
> (var(j).lt.param(2)) )
|
|
|
1467 |
> then
|
|
|
1468 |
istrue(j)=1
|
|
|
1469 |
else
|
|
|
1470 |
istrue(j)=0
|
|
|
1471 |
endif
|
|
|
1472 |
|
|
|
1473 |
c OUT
|
|
|
1474 |
elseif (icmd.eq.4) then
|
|
|
1475 |
if ( (var(j).lt.param(1)).or.
|
|
|
1476 |
> (var(j).gt.param(2)) )
|
|
|
1477 |
> then
|
|
|
1478 |
istrue(j)=1
|
|
|
1479 |
else
|
|
|
1480 |
istrue(j)=0
|
|
|
1481 |
endif
|
|
|
1482 |
|
|
|
1483 |
c EQ
|
|
|
1484 |
elseif (icmd.eq.5) then
|
|
|
1485 |
if (abs(var(j)-param(1)).lt.eps) then
|
|
|
1486 |
istrue(j)=1
|
|
|
1487 |
else
|
|
|
1488 |
istrue(j)=0
|
|
|
1489 |
endif
|
|
|
1490 |
|
|
|
1491 |
c TRUE
|
|
|
1492 |
elseif (icmd.eq.6) then
|
|
|
1493 |
if (abs(var(j)).lt.eps) then
|
|
|
1494 |
istrue(j)=0
|
|
|
1495 |
else
|
|
|
1496 |
istrue(j)=1
|
|
|
1497 |
endif
|
|
|
1498 |
|
|
|
1499 |
c FALSE
|
|
|
1500 |
elseif (icmd.eq.7) then
|
|
|
1501 |
if (abs(var(j)).lt.eps) then
|
|
|
1502 |
istrue(j)=1
|
|
|
1503 |
else
|
|
|
1504 |
istrue(j)=0
|
|
|
1505 |
endif
|
|
|
1506 |
|
|
|
1507 |
c ALL
|
|
|
1508 |
elseif (icmd.eq.8) then
|
|
|
1509 |
istrue(j) = 1
|
|
|
1510 |
do k=1,nval
|
|
|
1511 |
iparam = nint(param(k))-1
|
|
|
1512 |
if (btest(intvar(j),iparam).eqv..false.) then
|
|
|
1513 |
istrue(j) = 0
|
|
|
1514 |
endif
|
|
|
1515 |
enddo
|
|
|
1516 |
|
|
|
1517 |
c ANY
|
|
|
1518 |
elseif (icmd.eq.9) then
|
|
|
1519 |
istrue(j) = 0
|
|
|
1520 |
do k=1,nval
|
|
|
1521 |
iparam = nint(param(k))-1
|
|
|
1522 |
if (btest(intvar(j),iparam).eqv..true.) then
|
|
|
1523 |
istrue(j) = 1
|
|
|
1524 |
endif
|
|
|
1525 |
enddo
|
|
|
1526 |
|
|
|
1527 |
c NONE
|
|
|
1528 |
elseif (icmd.eq.10) then
|
|
|
1529 |
istrue(j) = 1
|
|
|
1530 |
do k=1,nval
|
|
|
1531 |
iparam = nint(param(k))-1
|
|
|
1532 |
if (btest(intvar(j),iparam).eqv..true.) then
|
|
|
1533 |
istrue(j) = 0
|
|
|
1534 |
endif
|
|
|
1535 |
enddo
|
|
|
1536 |
|
|
|
1537 |
endif
|
|
|
1538 |
|
|
|
1539 |
enddo
|
|
|
1540 |
|
|
|
1541 |
c --- Determine the overall boolean value ----------
|
|
|
1542 |
|
|
|
1543 |
c ALL
|
|
|
1544 |
if (itime.eq.1) then
|
|
|
1545 |
decision=1
|
|
|
1546 |
do j=1,nsel
|
|
|
1547 |
if (istrue(j).eq.0) then
|
|
|
1548 |
decision=0
|
|
|
1549 |
goto 110
|
|
|
1550 |
endif
|
|
|
1551 |
enddo
|
|
|
1552 |
110 continue
|
|
|
1553 |
|
|
|
1554 |
c ANY
|
|
|
1555 |
elseif (itime.eq.2) then
|
|
|
1556 |
decision=0
|
|
|
1557 |
do j=1,nsel
|
|
|
1558 |
if (istrue(j).eq.1) then
|
|
|
1559 |
decision=1
|
|
|
1560 |
goto 120
|
|
|
1561 |
endif
|
|
|
1562 |
enddo
|
|
|
1563 |
120 continue
|
|
|
1564 |
|
|
|
1565 |
c NONE
|
|
|
1566 |
elseif (itime.eq.3) then
|
|
|
1567 |
decision=1
|
|
|
1568 |
do j=1,nsel
|
|
|
1569 |
if (istrue(j).eq.1) then
|
|
|
1570 |
decision=0
|
|
|
1571 |
goto 130
|
|
|
1572 |
endif
|
|
|
1573 |
enddo
|
|
|
1574 |
130 continue
|
|
|
1575 |
|
|
|
1576 |
c TRIGGER
|
|
|
1577 |
elseif (itime.lt.0) then
|
|
|
1578 |
decision=1
|
|
|
1579 |
do j=1,nsel
|
|
|
1580 |
if (istrue(j).eq.1) then
|
|
|
1581 |
trigger(j) = ior( trigger(j), 2**(abs(itime)-1) )
|
|
|
1582 |
endif
|
|
|
1583 |
enddo
|
|
|
1584 |
|
|
|
1585 |
endif
|
|
|
1586 |
|
|
|
1587 |
c --- Put the new boolean value onto the stack
|
|
|
1588 |
|
|
|
1589 |
nstack=nstack+1
|
|
|
1590 |
stack(nstack)=decision
|
|
|
1591 |
|
|
|
1592 |
c Exit point for loop
|
|
|
1593 |
200 continue
|
|
|
1594 |
goto 100
|
|
|
1595 |
|
|
|
1596 |
endif
|
|
|
1597 |
|
|
|
1598 |
c Return the decision (selected or non-selected)
|
|
|
1599 |
300 continue
|
|
|
1600 |
|
|
|
1601 |
select=stack(1)
|
|
|
1602 |
|
|
|
1603 |
end
|
|
|
1604 |
|
|
|
1605 |
|
|
|
1606 |
c --------------------------------------------------------------------------
|
|
|
1607 |
c Split a region string and get corners of the domain
|
|
|
1608 |
c --------------------------------------------------------------------------
|
|
|
1609 |
|
|
|
1610 |
subroutine regionsplit(string,iregion,xcorner,ycorner)
|
|
|
1611 |
|
|
|
1612 |
c The region string comes either as <lonw,lone,lats,latn> or as <lon1,lat1,
|
|
|
1613 |
c lon2,lat2,lon3,lat3,lon4,lat4>: split it into ints components and get the
|
|
|
1614 |
c four coordinates for the region
|
|
|
1615 |
|
|
|
1616 |
implicit none
|
|
|
1617 |
|
|
|
1618 |
c Declaration of subroutine parameters
|
|
|
1619 |
character*80 string
|
|
|
1620 |
real xcorner(4),ycorner(4)
|
|
|
1621 |
integer iregion
|
|
|
1622 |
|
|
|
1623 |
c Local variables
|
|
|
1624 |
integer i,n
|
|
|
1625 |
integer il,ir
|
|
|
1626 |
real subfloat (80)
|
|
|
1627 |
integer stat
|
|
|
1628 |
integer len
|
|
|
1629 |
|
|
|
1630 |
c ------- Split the string
|
|
|
1631 |
i = 1
|
|
|
1632 |
n = 0
|
|
|
1633 |
stat = 0
|
|
|
1634 |
il = 1
|
|
|
1635 |
len = len_trim(string)
|
|
|
1636 |
|
|
|
1637 |
100 continue
|
|
|
1638 |
|
|
|
1639 |
c Find start of a substring
|
|
|
1640 |
do while ( stat.eq.0 )
|
|
|
1641 |
if ( string(i:i).ne.' ' ) then
|
|
|
1642 |
stat = 1
|
|
|
1643 |
il = i
|
|
|
1644 |
else
|
|
|
1645 |
i = i + 1
|
|
|
1646 |
endif
|
|
|
1647 |
enddo
|
|
|
1648 |
|
|
|
1649 |
c Find end of substring
|
|
|
1650 |
do while ( stat.eq.1 )
|
|
|
1651 |
if ( ( string(i:i).eq.' ' ) .or. ( i.eq.len ) ) then
|
|
|
1652 |
stat = 2
|
|
|
1653 |
ir = i
|
|
|
1654 |
else
|
|
|
1655 |
i = i + 1
|
|
|
1656 |
endif
|
|
|
1657 |
enddo
|
|
|
1658 |
|
|
|
1659 |
c Convert the substring into a number
|
|
|
1660 |
if ( stat.eq.2 ) then
|
|
|
1661 |
n = n + 1
|
|
|
1662 |
read(string(il:ir),*) subfloat(n)
|
|
|
1663 |
stat = 0
|
|
|
1664 |
endif
|
|
|
1665 |
|
|
|
1666 |
if ( i.lt.len ) goto 100
|
|
|
1667 |
|
|
|
1668 |
|
|
|
1669 |
c -------- Get the region number
|
|
|
1670 |
|
|
|
1671 |
iregion = nint(subfloat(1))
|
|
|
1672 |
|
|
|
1673 |
c -------- Get the corners of the region
|
|
|
1674 |
|
|
|
1675 |
if ( n.eq.5 ) then ! lonw(2),lone(3),lats(4),latn(5)
|
|
|
1676 |
|
|
|
1677 |
xcorner(1) = subfloat(2)
|
|
|
1678 |
ycorner(1) = subfloat(4)
|
|
|
1679 |
|
|
|
1680 |
xcorner(2) = subfloat(3)
|
|
|
1681 |
ycorner(2) = subfloat(4)
|
|
|
1682 |
|
|
|
1683 |
xcorner(3) = subfloat(3)
|
|
|
1684 |
ycorner(3) = subfloat(5)
|
|
|
1685 |
|
|
|
1686 |
xcorner(4) = subfloat(2)
|
|
|
1687 |
ycorner(4) = subfloat(5)
|
|
|
1688 |
|
|
|
1689 |
elseif ( n.eq.9 ) then ! lon1,lat1,lon2,lat2,lon3,lon4,lat4
|
|
|
1690 |
|
|
|
1691 |
xcorner(1) = subfloat(2)
|
|
|
1692 |
ycorner(1) = subfloat(3)
|
|
|
1693 |
|
|
|
1694 |
xcorner(2) = subfloat(4)
|
|
|
1695 |
ycorner(2) = subfloat(5)
|
|
|
1696 |
|
|
|
1697 |
xcorner(3) = subfloat(6)
|
|
|
1698 |
ycorner(3) = subfloat(7)
|
|
|
1699 |
|
|
|
1700 |
xcorner(4) = subfloat(8)
|
|
|
1701 |
ycorner(4) = subfloat(9)
|
|
|
1702 |
|
|
|
1703 |
else
|
|
|
1704 |
|
|
|
1705 |
print*,' ERROR: invalid region specification '
|
|
|
1706 |
print*,' ',trim(string)
|
|
|
1707 |
stop
|
|
|
1708 |
|
|
|
1709 |
endif
|
|
|
1710 |
|
|
|
1711 |
|
|
|
1712 |
end
|
|
|
1713 |
|
|
|
1714 |
c --------------------------------------------------------------------------
|
|
|
1715 |
c Decide whether lat/lon point is in or out of region
|
|
|
1716 |
c --------------------------------------------------------------------------
|
|
|
1717 |
|
|
|
1718 |
integer function inregion (lon,lat,xcorner,ycorner)
|
|
|
1719 |
|
|
|
1720 |
c Decide whether point (lon/lat) is in the region specified by <xcorner(1..4),
|
|
|
1721 |
c ycorner(1..4).
|
|
|
1722 |
|
|
|
1723 |
implicit none
|
|
|
1724 |
|
|
|
1725 |
c Declaration of subroutine parameters
|
|
|
1726 |
real lon,lat
|
|
|
1727 |
real xcorner(4),ycorner(4)
|
|
|
1728 |
|
|
|
1729 |
c Local variables
|
|
|
1730 |
integer flag
|
|
|
1731 |
real xmin,xmax,ymin,ymax
|
|
|
1732 |
integer i
|
|
|
1733 |
|
|
|
1734 |
c Reset the flag
|
|
|
1735 |
flag = 0
|
|
|
1736 |
|
|
|
1737 |
c Set some boundaries
|
|
|
1738 |
xmax = xcorner(1)
|
|
|
1739 |
xmin = xcorner(1)
|
|
|
1740 |
ymax = ycorner(1)
|
|
|
1741 |
ymin = ycorner(1)
|
|
|
1742 |
do i=2,4
|
|
|
1743 |
if (xcorner(i).lt.xmin) xmin = xcorner(i)
|
|
|
1744 |
if (xcorner(i).gt.xmax) xmax = xcorner(i)
|
|
|
1745 |
if (ycorner(i).lt.ymin) ymin = ycorner(i)
|
|
|
1746 |
if (ycorner(i).gt.ymax) ymax = ycorner(i)
|
|
|
1747 |
enddo
|
|
|
1748 |
|
|
|
1749 |
c Do the tests - set flag=1 if all tests pased
|
|
|
1750 |
if (lon.lt.xmin) goto 970
|
|
|
1751 |
if (lon.gt.xmax) goto 970
|
|
|
1752 |
if (lat.lt.ymin) goto 970
|
|
|
1753 |
if (lat.gt.ymax) goto 970
|
|
|
1754 |
|
|
|
1755 |
if ((lon-xcorner(1))*(ycorner(2)-ycorner(1))-
|
|
|
1756 |
> (lat-ycorner(1))*(xcorner(2)-xcorner(1)).gt.0.) goto 970
|
|
|
1757 |
if ((lon-xcorner(2))*(ycorner(3)-ycorner(2))-
|
|
|
1758 |
> (lat-ycorner(2))*(xcorner(3)-xcorner(2)).gt.0.) goto 970
|
|
|
1759 |
if ((lon-xcorner(3))*(ycorner(4)-ycorner(3))-
|
|
|
1760 |
> (lat-ycorner(3))*(xcorner(4)-xcorner(3)).gt.0.) goto 970
|
|
|
1761 |
if ((lon-xcorner(4))*(ycorner(1)-ycorner(4))-
|
|
|
1762 |
> (lat-ycorner(4))*(xcorner(1)-xcorner(4)).gt.0.) goto 970
|
|
|
1763 |
|
|
|
1764 |
flag = 1
|
|
|
1765 |
|
|
|
1766 |
c Return the value
|
|
|
1767 |
970 continue
|
|
|
1768 |
|
|
|
1769 |
inregion = flag
|
|
|
1770 |
|
|
|
1771 |
return
|
|
|
1772 |
|
|
|
1773 |
end
|
|
|
1774 |
|
|
|
1775 |
|
|
|
1776 |
c --------------------------------------------------------------------------
|
|
|
1777 |
c Spherical distance between lat/lon points
|
|
|
1778 |
c --------------------------------------------------------------------------
|
|
|
1779 |
|
|
|
1780 |
real function sdis(xp,yp,xq,yq)
|
|
|
1781 |
c
|
|
|
1782 |
c calculates spherical distance (in km) between two points given
|
|
|
1783 |
c by their spherical coordinates (xp,yp) and (xq,yq), respectively.
|
|
|
1784 |
c
|
|
|
1785 |
real re
|
|
|
1786 |
parameter (re=6370.)
|
|
|
1787 |
real pi180
|
|
|
1788 |
parameter (pi180=3.14159/180.)
|
|
|
1789 |
real xp,yp,xq,yq,arg
|
|
|
1790 |
|
|
|
1791 |
arg=sin(pi180*yp)*sin(pi180*yq)+
|
|
|
1792 |
> cos(pi180*yp)*cos(pi180*yq)*cos(pi180*(xp-xq))
|
|
|
1793 |
if (arg.lt.-1.) arg=-1.
|
|
|
1794 |
if (arg.gt.1.) arg=1.
|
|
|
1795 |
|
|
|
1796 |
sdis=re*acos(arg)
|
|
|
1797 |
|
|
|
1798 |
end
|
|
|
1799 |
|
|
|
1800 |
|
|
|
1801 |
c ****************************************************************
|
|
|
1802 |
c * Given some spherical polygon S and some point X known to be *
|
|
|
1803 |
c * located inside S, these routines will determine if an arbit- *
|
|
|
1804 |
c * -rary point P lies inside S, outside S, or on its boundary. *
|
|
|
1805 |
c * The calling program must first call DefSPolyBndry to define *
|
|
|
1806 |
c * the boundary of S and the point X. Any subsequent call to *
|
|
|
1807 |
c * subroutine LctPtRelBndry will determine if some point P lies *
|
|
|
1808 |
c * inside or outside S, or on its boundary. (Usually *
|
|
|
1809 |
c * DefSPolyBndry is called once, then LctPrRelBndry is called *
|
|
|
1810 |
c * many times). *
|
|
|
1811 |
c * *
|
|
|
1812 |
c * REFERENCE: Bevis, M. and Chatelain, J.-L. (1989) *
|
|
|
1813 |
c * Maflaematical Geology, vol 21. *
|
|
|
1814 |
c * VERSION 1.0 *
|
|
|
1815 |
c ****************************************************************
|
|
|
1816 |
|
|
|
1817 |
Subroutine DefSPolyBndry(vlat,vlon,nv,xlat, xlon)
|
|
|
1818 |
|
|
|
1819 |
c ****************************************************************
|
|
|
1820 |
c * This mmn entry point is used m define ~e spheric~ polygon S *
|
|
|
1821 |
c * and the point X. *
|
|
|
1822 |
c * ARGUMENTS: *
|
|
|
1823 |
c * vlat,vlon (sent) ... vectors containing the latitude and *
|
|
|
1824 |
c * longitude of each vertex of the *
|
|
|
1825 |
c * spherical polygon S. The ith.vertex is *
|
|
|
1826 |
c * located at [vlat(i),vlon(i)]. *
|
|
|
1827 |
c * nv (sent) ... the number of vertices and sides in the *
|
|
|
1828 |
c * spherical polygon S *
|
|
|
1829 |
c * xlat,xlon (sent) ... latitude and longitude of some point X *
|
|
|
1830 |
c * located inside S. X must not be located *
|
|
|
1831 |
c * on any great circle that includes two *
|
|
|
1832 |
c * vertices of S. *
|
|
|
1833 |
c * *
|
|
|
1834 |
c * UNITS AND SIGN CONVENTION: *
|
|
|
1835 |
c * Latitudes and longitudes are specified in degrees. *
|
|
|
1836 |
c * Latitudes are positive to the north and negative to the *
|
|
|
1837 |
c * south. *
|
|
|
1838 |
c * Longitudes are positive to the east and negative to the *
|
|
|
1839 |
c * west. *
|
|
|
1840 |
c * *
|
|
|
1841 |
c * VERTEX ENUMERATION: *
|
|
|
1842 |
c * The vertices of S should be numbered sequentially around the *
|
|
|
1843 |
c * border of the spherical polygon. Vertex 1 lies between vertex*
|
|
|
1844 |
c * nv and vertex 2. Neighbouring vertices must be seperated by *
|
|
|
1845 |
c * less than 180 degrees. (In order to generate a polygon side *
|
|
|
1846 |
c * whose arc length equals or exceeds 180 degrees simply *
|
|
|
1847 |
c * introduce an additional (pseudo)vertex). Having chosen *
|
|
|
1848 |
c * vertex 1, the user may number the remaining vertices in *
|
|
|
1849 |
c * either direction. However if the user wishes to use the *
|
|
|
1850 |
c * subroutine SPA to determine the area of the polygon S (Bevis *
|
|
|
1851 |
c * & Cambareri, 1987, Math. Geol., v.19, p. 335-346) then he or *
|
|
|
1852 |
c * she must follow the convention whereby in moving around the *
|
|
|
1853 |
c * polygon border in the direction of increasing vertex number *
|
|
|
1854 |
c * clockwise bends occur at salient vertices. A vertex is *
|
|
|
1855 |
c * salient if the interior angle is less than 180 degrees. *
|
|
|
1856 |
c * (In the case of a convex polygon this convention implies *
|
|
|
1857 |
c * that vertices are numbered in clockwise sequence). *
|
|
|
1858 |
c ****************************************************************
|
|
|
1859 |
|
|
|
1860 |
implicit none
|
|
|
1861 |
|
|
|
1862 |
integer mxnv,nv
|
|
|
1863 |
|
|
|
1864 |
c ----------------------------------------------------------------
|
|
|
1865 |
c Edit next statement to increase maximum number of vertices that
|
|
|
1866 |
c may be used to define the spherical polygon S
|
|
|
1867 |
c The value of parameter mxnv in subroutine LctPtRelBndry must match
|
|
|
1868 |
c that of parameter mxnv in this subroutine, as assigned above.
|
|
|
1869 |
c ----------------------------------------------------------------
|
|
|
1870 |
parameter (mxnv=2000)
|
|
|
1871 |
|
|
|
1872 |
real vlat(nv),vlon(nv),xlat,xlon,dellon
|
|
|
1873 |
real tlonv(mxnv),vlat_c(mxnv),vlon_c(mxnv),xlat_c,xlon_c
|
|
|
1874 |
integer i,ibndry,nv_c,ip
|
|
|
1875 |
|
|
|
1876 |
data ibndry/0/
|
|
|
1877 |
|
|
|
1878 |
common/spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry
|
|
|
1879 |
|
|
|
1880 |
if (nv.gt.mxnv) then
|
|
|
1881 |
print *,'nv exceeds maximum allowed value'
|
|
|
1882 |
print *,'adjust parameter mxnv in subroutine DefSPolyBndry'
|
|
|
1883 |
stop
|
|
|
1884 |
endif
|
|
|
1885 |
|
|
|
1886 |
ibndry=1 ! boundary defined at least once (flag)
|
|
|
1887 |
nv_c=nv ! copy for named common
|
|
|
1888 |
xlat_c=xlat ! . . . .
|
|
|
1889 |
xlon_c=xlon !
|
|
|
1890 |
|
|
|
1891 |
do i=1,nv
|
|
|
1892 |
vlat_c(i)=vlat(i) ! "
|
|
|
1893 |
vlon_c(i)=vlon(i) !
|
|
|
1894 |
|
|
|
1895 |
call TrnsfmLon(xlat,xlon,vlat(i),vlon(i),tlonv(i))
|
|
|
1896 |
|
|
|
1897 |
if (i.gt.1) then
|
|
|
1898 |
ip=i-1
|
|
|
1899 |
else
|
|
|
1900 |
ip=nv
|
|
|
1901 |
endif
|
|
|
1902 |
|
|
|
1903 |
if ((vlat(i).eq.vlat(ip)).and.(vlon(i).eq.vlon(ip))) then
|
|
|
1904 |
print *,'DefSPolyBndry detects user error:'
|
|
|
1905 |
print *,'vertices ',i,' and ',ip,' are not distinct'
|
|
|
1906 |
print*,'lat ',i,ip,vlat(i),vlat(ip)
|
|
|
1907 |
print*,'lon ',i,ip,vlon(i),vlon(ip)
|
|
|
1908 |
stop
|
|
|
1909 |
endif
|
|
|
1910 |
|
|
|
1911 |
if (tlonv(i).eq.tlonv(ip)) then
|
|
|
1912 |
print *,'DefSPolyBndry detects user error:'
|
|
|
1913 |
print *,'vertices ',i,' & ',ip,' on same gt. circle as X'
|
|
|
1914 |
stop
|
|
|
1915 |
endif
|
|
|
1916 |
|
|
|
1917 |
if (vlat(i).eq.(-vlat(ip))) then
|
|
|
1918 |
dellon=vlon(i)-vlon(ip)
|
|
|
1919 |
if (dellon.gt.+180.) dellon=dellon-360.
|
|
|
1920 |
if (dellon.lt.-180.) dellon=dellon-360.
|
|
|
1921 |
if ((dellon.eq.+180.0).or.(dellon.eq.-180.0)) then
|
|
|
1922 |
print *,'DefSPolyBndry detects user error:'
|
|
|
1923 |
print *,'vertices ',i,' and ',ip,' are antipodal'
|
|
|
1924 |
stop
|
|
|
1925 |
endif
|
|
|
1926 |
endif
|
|
|
1927 |
enddo
|
|
|
1928 |
|
|
|
1929 |
return
|
|
|
1930 |
|
|
|
1931 |
end
|
|
|
1932 |
|
|
|
1933 |
|
|
|
1934 |
c ****************************************************************
|
|
|
1935 |
|
|
|
1936 |
Subroutine LctPtRelBndry(plat,plon,location)
|
|
|
1937 |
|
|
|
1938 |
c ****************************************************************
|
|
|
1939 |
|
|
|
1940 |
c ****************************************************************
|
|
|
1941 |
c * This routine is used to see if some point P is located *
|
|
|
1942 |
c * inside, outside or on the boundary of the spherical polygon *
|
|
|
1943 |
c * S previously defined by a call to subroutine DefSPolyBndry. *
|
|
|
1944 |
c * There is a single restriction on point P: it must not be *
|
|
|
1945 |
c * antipodal to the point X defined in the call to DefSPolyBndry*
|
|
|
1946 |
c * (ie.P and X cannot be seperated by exactly 180 degrees). *
|
|
|
1947 |
c * ARGUMENTS: *
|
|
|
1948 |
c * plat,plon (sent)... the latitude and longitude of point P *
|
|
|
1949 |
c * location (returned)... specifies the location of P: *
|
|
|
1950 |
c * location=0 implies P is outside of S *
|
|
|
1951 |
c * location=1 implies P is inside of S *
|
|
|
1952 |
c * location=2 implies P on boundary of S *
|
|
|
1953 |
c * location=3 implies user error (P is *
|
|
|
1954 |
c * antipodal to X) *
|
|
|
1955 |
c * UNFfS AND SIGN CONVENTION: *
|
|
|
1956 |
c * Latitudes and longitudes are specified in degrees. *
|
|
|
1957 |
c * Latitudes are positive to the north and negative to the *
|
|
|
1958 |
c * south. *
|
|
|
1959 |
c * Longitudes are positive to the east and negative to the *
|
|
|
1960 |
c * west. *
|
|
|
1961 |
c ****************************************************************
|
|
|
1962 |
|
|
|
1963 |
implicit none
|
|
|
1964 |
|
|
|
1965 |
integer mxnv
|
|
|
1966 |
|
|
|
1967 |
c ----------------------------------------------------------------
|
|
|
1968 |
c The statement below must match that in subroutine DefSPolyBndry
|
|
|
1969 |
c ----------------------------------------------------------------
|
|
|
1970 |
|
|
|
1971 |
parameter (mxnv=2000)
|
|
|
1972 |
|
|
|
1973 |
real tlonv(mxnv),vlat_c(mxnv),vlon_c(mxnv),xlat_c,xlon_c
|
|
|
1974 |
real plat,plon,vAlat,vAlon,vBlat,vBlon,tlonA,tlonB,tlonP
|
|
|
1975 |
real tlon_X,tlon_P,tlon_B,dellon
|
|
|
1976 |
integer i,ibndry,nv_c,location,icross,ibrngAB,ibrngAP,ibrngPB
|
|
|
1977 |
integer ibrng_BX,ibrng_BP,istrike
|
|
|
1978 |
|
|
|
1979 |
common/spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry
|
|
|
1980 |
|
|
|
1981 |
if (ibndry.eq.0) then ! user has never defined the bndry
|
|
|
1982 |
print*,'Subroutine LctPtRelBndry detects user error:'
|
|
|
1983 |
print*,'Subroutine DefSPolyBndry must be called before'
|
|
|
1984 |
print*,'subroutine LctPtRelBndry can be called'
|
|
|
1985 |
stop
|
|
|
1986 |
endif
|
|
|
1987 |
|
|
|
1988 |
if (plat.eq.(-xlat_c)) then
|
|
|
1989 |
dellon=plon-xlon_c
|
|
|
1990 |
if (dellon.lt.(-180.)) dellon=dellon+360.
|
|
|
1991 |
if (dellon.gt.+180.) dellon=dellon-360.
|
|
|
1992 |
if ((dellon.eq.+180.0).or.(dellon.eq.-180.)) then
|
|
|
1993 |
print*,'Warning: LctPtRelBndry detects case P antipodal
|
|
|
1994 |
> to X'
|
|
|
1995 |
print*,'location of P relative to S is undetermined'
|
|
|
1996 |
location=3
|
|
|
1997 |
return
|
|
|
1998 |
endif
|
|
|
1999 |
endif
|
|
|
2000 |
|
|
|
2001 |
location=0 ! default ( P is outside S)
|
|
|
2002 |
icross=0 ! initialize counter
|
|
|
2003 |
|
|
|
2004 |
if ((plat.eq.xlat_c).and.(plon.eq.xlon_c)) then
|
|
|
2005 |
location=1
|
|
|
2006 |
return
|
|
|
2007 |
endif
|
|
|
2008 |
|
|
|
2009 |
|
|
|
2010 |
call TrnsfmLon (xlat_c,xlon_c,plat,plon,tlonP)
|
|
|
2011 |
|
|
|
2012 |
do i=1,nv_c ! start of loop over sides of S
|
|
|
2013 |
|
|
|
2014 |
vAlat=vlat_c(i)
|
|
|
2015 |
vAlon=vlon_c(i)
|
|
|
2016 |
tlonA=tlonv(i)
|
|
|
2017 |
|
|
|
2018 |
if (i.lt.nv_c) then
|
|
|
2019 |
vBlat=vlat_c(i+1)
|
|
|
2020 |
vBlon=vlon_c(i+1)
|
|
|
2021 |
tlonB=tlonv(i+1)
|
|
|
2022 |
else
|
|
|
2023 |
vBlat=vlat_c(1)
|
|
|
2024 |
vBlon=vlon_c(1)
|
|
|
2025 |
tlonB=tlonv(1)
|
|
|
2026 |
endif
|
|
|
2027 |
|
|
|
2028 |
istrike=0
|
|
|
2029 |
|
|
|
2030 |
if (tlonP.eq.tlonA) then
|
|
|
2031 |
istrike=1
|
|
|
2032 |
else
|
|
|
2033 |
call EastOrWest(tlonA,tlonB,ibrngAB)
|
|
|
2034 |
call EastOrWest(tlonA,tlonP,ibrngAP)
|
|
|
2035 |
call EastOrWest(tlonP,tlonB,ibrngPB)
|
|
|
2036 |
|
|
|
2037 |
|
|
|
2038 |
if((ibrngAP.eq.ibrngAB).and.(ibrngPB.eq.ibrngAB)) istrike=1
|
|
|
2039 |
endif
|
|
|
2040 |
|
|
|
2041 |
|
|
|
2042 |
if (istrike.eq.1) then
|
|
|
2043 |
|
|
|
2044 |
if ((plat.eq.vAlat).and.(plon.eq.vAlon)) then
|
|
|
2045 |
location=2 ! P lies on a vertex of S
|
|
|
2046 |
return
|
|
|
2047 |
endif
|
|
|
2048 |
call TrnsfmLon(vAlat,vAlon,xlat_c,xlon_c,tlon_X)
|
|
|
2049 |
call TrnsfmLon(vAlat,vAlon,vBlat,vBlon,tlon_B)
|
|
|
2050 |
call TrnsfmLon(vAlat,vAlon,plat,plon,tlon_P)
|
|
|
2051 |
|
|
|
2052 |
if (tlon_P.eq.tlon_B) then
|
|
|
2053 |
location=2 ! P lies on side of S
|
|
|
2054 |
return
|
|
|
2055 |
else
|
|
|
2056 |
call EastOrWest(tlon_B,tlon_X,ibrng_BX)
|
|
|
2057 |
call EastOrWest(tlon_B,tlon_P,ibrng_BP)
|
|
|
2058 |
if(ibrng_BX.eq.(-ibrng_BP)) icross=icross+1
|
|
|
2059 |
endif
|
|
|
2060 |
|
|
|
2061 |
endif
|
|
|
2062 |
enddo ! end of loop over the sides of S
|
|
|
2063 |
|
|
|
2064 |
|
|
|
2065 |
c if the arc XP crosses the boundary S an even number of times then P
|
|
|
2066 |
c is in S
|
|
|
2067 |
|
|
|
2068 |
if (mod(icross,2).eq.0) location=1
|
|
|
2069 |
|
|
|
2070 |
return
|
|
|
2071 |
|
|
|
2072 |
end
|
|
|
2073 |
|
|
|
2074 |
|
|
|
2075 |
c ****************************************************************
|
|
|
2076 |
|
|
|
2077 |
subroutine TrnsfmLon(plat,plon,qlat,qlon,tranlon)
|
|
|
2078 |
|
|
|
2079 |
c ****************************************************************
|
|
|
2080 |
c * This subroutine is required by subroutines DefSPolyBndry & *
|
|
|
2081 |
c * LctPtRelBndry. It finds the 'longitude' of point Q in a *
|
|
|
2082 |
c * geographic coordinate system for which point P acts as a *
|
|
|
2083 |
c * 'north pole'. SENT: plat,plon,qlat,qlon, in degrees. *
|
|
|
2084 |
c * RETURNED: tranlon, in degrees. *
|
|
|
2085 |
c ****************************************************************
|
|
|
2086 |
|
|
|
2087 |
implicit none
|
|
|
2088 |
|
|
|
2089 |
real pi,dtr,plat,plon,qlat,qlon,tranlon,t,b
|
|
|
2090 |
parameter (pi=3.141592654,dtr=pi/180.0)
|
|
|
2091 |
|
|
|
2092 |
if (plat.eq.90.) then
|
|
|
2093 |
tranlon=qlon
|
|
|
2094 |
else
|
|
|
2095 |
t=sin((qlon-plon)*dtr)*cos(qlat*dtr)
|
|
|
2096 |
b=sin(dtr*qlat)*cos(plat*dtr)-cos(qlat*dtr)*sin(plat*dtr)
|
|
|
2097 |
> *cos((qlon-plon)*dtr)
|
|
|
2098 |
tranlon=atan2(t,b)/dtr
|
|
|
2099 |
endif
|
|
|
2100 |
|
|
|
2101 |
return
|
|
|
2102 |
end
|
|
|
2103 |
|
|
|
2104 |
c ****************************************************************
|
|
|
2105 |
|
|
|
2106 |
subroutine EastOrWest(clon,dlon,ibrng)
|
|
|
2107 |
|
|
|
2108 |
c ****************************************************************
|
|
|
2109 |
c * This subroutine is required by subroutine LctPtRelBndry. *
|
|
|
2110 |
c * This routine determines if in travelling the shortest path *
|
|
|
2111 |
c * from point C (at longitude clon) to point D (at longitude *
|
|
|
2112 |
c * dlon) one is heading east, west or neither. *
|
|
|
2113 |
c * SENT: clon,dlon; in degrees. RETURNED: ibrng *
|
|
|
2114 |
c * (1=east,-1=west, 0=neither). *
|
|
|
2115 |
c ****************************************************************
|
|
|
2116 |
|
|
|
2117 |
implicit none
|
|
|
2118 |
real clon,dlon,del
|
|
|
2119 |
integer ibrng
|
|
|
2120 |
del=dlon-clon
|
|
|
2121 |
if (del.gt.180.) del=del-360.
|
|
|
2122 |
if (del.lt.-180.) del=del+360.
|
|
|
2123 |
if ((del.gt.0.0).and.(del.ne.180.)) then
|
|
|
2124 |
ibrng=-1 ! (D is west of C)
|
|
|
2125 |
elseif ((del.lt.0.0).and.(del.ne.-180.)) then
|
|
|
2126 |
ibrng=+1 ! (D is east of C)
|
|
|
2127 |
else
|
|
|
2128 |
ibrng=0 ! (D north or south of C)
|
|
|
2129 |
endif
|
|
|
2130 |
return
|
|
|
2131 |
end
|