Subversion Repositories lagranto.20cr

Rev

Details | Last modification | View Log | RSS feed

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