Subversion Repositories lagranto.ecmwf

Rev

Rev 15 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

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