Subversion Repositories lagranto.arpege

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 michaesp 1
      PROGRAM create_startf
2
 
3
c     **************************************************************
4
c     * Create a <startfile> for <lagrangto>. It can be chosen     *
5
c     * whether to start from an isentropic or an isobaric         *
6
c     * surface. The starting points are equidistantly distributed *
7
c     * Michael Sprenger / Autumn 2004                             *
8
c     **************************************************************
9
 
10
      use netcdf
11
 
12
      implicit none
13
 
14
c     --------------------------------------------------------------
15
c     Set parameters
16
c     --------------------------------------------------------------
17
 
18
c     Maximum number of starting positions
19
      integer           nmax
20
      parameter         (nmax=10000000)
21
 
22
c     Maximum number of model levels
23
      integer           nlevmax
24
      parameter         (nlevmax=200)
25
 
26
c     Grid constant (distance in km corresponding to 1 deg at the equator)
27
      real              deltat
28
      parameter         (deltat=111.)
29
 
30
c     Mathematical constant (conversion degree -> radian)
31
      real              pi180
32
      parameter         (pi180=3.14159/180.)
33
 
34
c     Numerical epsilon
35
      real              eps
36
      parameter         (eps=0.00001)
37
 
38
c     --------------------------------------------------------------
39
c     Set variables
40
c     --------------------------------------------------------------
41
 
42
c     Filenames and output format
43
      character*80      pfile0,pfile1                  ! P filenames
44
      character*80      sfile0,sfile1                  ! S filenames
45
      character*80      ofile                          ! Output filename
46
      integer           oformat                        ! Output format
47
      real              timeshift                      ! Time shift relative to data files <*0>
48
      real              timeinc                        ! Time increment between input files
49
 
50
c     Horizontal grid
51
      character*80      hmode                          ! Horizontale mode
52
      real              lat1,lat2,lon1,lon2            ! Lat/lon boundaries
53
      real              ds,dlon,dlat                   ! Distance and lat/lon shifts
54
      character*80      hfile                          ! Filename
55
      integer           hn                             ! Number of entries in lat/lon list 
56
      real              latlist(nmax)                  ! List of latitudes
57
      real              lonlist(nmax)                  ! List of longitudes
58
      integer           pn                             ! Number of entries in lat/lon poly
59
      real              latpoly(500)                   ! List of polygon latitudes
60
      real              lonpoly(500)                   ! List of polygon longitudes
61
      real              loninpoly,latinpoly            ! Lon/lat inside polygon
62
      character*80      regionf                        ! Region file
63
      integer           iregion                        ! Region number
64
      real              xcorner(4),ycorner(4)          ! Vertices of region
65
 
66
c     Vertical grid
67
      character*80      vmode                          ! Vertical mode
68
      real              lev1,lev2,levlist(nmax)        ! Single levels, and list of levels
69
      character*80      vfile                          ! Filename
70
      integer           vn                             ! Number of entries
71
 
72
c     Unit of vertical axis
73
      character*80      umode                          ! Unit of vertical axis
74
 
75
c     Flag for 'no time check'
76
      character*80      timecheck                      ! Either 'no' or 'yes'
77
 
78
c     List of all starting positions
79
      integer           start_n                        ! Number of coordinates
80
      real              start_lat(nmax)                ! Latitudes
81
      real              start_lon(nmax)                ! Longitudes
82
      real              start_lev(nmax)                ! Levels (depending on vertical unit)
83
      real              start_pre(nmax)                ! Level in hPa
84
      integer           start_valid(nmax)              ! Flag whether starting height is valid
85
      integer           reftime(6)                     ! Reference time
86
      character*80      vars(10)                       ! Name of output fields (time,lon,lat,p)
87
      real,allocatable, dimension (:,:,:) :: tra       ! Trajectories (ntra,ntim,ncol)
88
      real              latmin,latmax
89
      real              lonmin,lonmax
90
      real              premin,premax
91
 
92
c     Grid description
93
      real              pollon,pollat                  ! Longitude/latitude of pole
94
      real              ak(nlevmax)                    ! Vertical layers and levels
95
      real              bk(nlevmax)
96
      real              xmin,xmax                      ! Zonal grid extension
97
      real              ymin,ymax                      ! Meridional grid extension
98
      integer           nx,ny,nz                       ! Grid dimensions
99
      real              dx,dy                          ! Horizontal grid resolution
100
      real,allocatable, dimension (:,:,:) :: pr        ! 3d pressure
101
      real,allocatable, dimension (:,:)   :: prs       ! surface pressure
102
      real,allocatable, dimension (:,:,:) :: th        ! 3d potential temperature
103
      real,allocatable, dimension (:,:)   :: ths       ! surface poential temperature
104
      real,allocatable, dimension (:,:,:) :: pv        ! 3d potential vorticity
105
      real,allocatable, dimension (:,:)   :: pvs       ! surface potential vorticiy
106
      real,allocatable, dimension (:,:,:) :: in        ! 3d 'dummy' array with vertical indices
107
      character*80      varname                        ! Name of input variable
108
      integer           fid                            ! File identifier
109
      real              stagz                          ! Vertical staggering
110
      real              mdv                            ! Missing data values
111
      real              tstart,tend                    ! Time on P and S file
112
      real              rid,rjd,rkd                    ! Real grid position
113
 
114
c     Auxiliary variable
115
      integer           i,j,k
116
      real              lon,lat
117
      real              rd
118
      integer           stat,flag
119
      real              tmp1,tmp2
120
      real              tfrac,frac
121
      real              radius,dist
122
      character*80      string
123
      character*80      selectstr
124
      character*80      umode_save
125
      character*80      maskname
126
      real              maskvalue
127
      character*2       maskoper
128
      integer           cdfid,varid,ierr
129
      integer           nmask
130
      integer           indx,indy
131
      character         ch1,ch2
132
      integer           len
133
      real,allocatable, dimension (:,:,:) :: fld0
134
      real,allocatable, dimension (:,:,:) :: fld1
135
      real,allocatable, dimension (:,:  ) :: sfc0
136
      real,allocatable, dimension (:,:)   :: sfc1
137
      real,allocatable, dimension (:,:)   :: mask
138
      integer           n_outside
139
 
140
c     Externals 
141
      real              int_index3     ! 3d interpolation
142
      external          int_index3
143
      real              sdis           ! Speherical distance
144
      external          sdis
145
      integer           inregion       ! In/out of region
146
      external          inrehion
147
 
148
c     ------------------------------------------------------------------
149
c     Start of program, Read parameters
150
c     ------------------------------------------------------------------
151
 
152
c     Write start message
153
      print*,'========================================================='
154
      print*,'         *** START OF PROGRAM CREATE_STARTF ***'
155
      print*
156
 
157
c     Read parameter file
158
      open(10,file='create_startf.param')
159
 
160
c      Input P and S file
161
       read(10,*) pfile0,pfile1
162
       read(10,*) sfile0,sfile1
163
       read(10,*) ofile
164
 
165
c      Read name of region file
166
       read(10,*) regionf
167
 
168
c      Reference time
169
       do i=1,6
170
          read(10,*) reftime(i)
171
       enddo
172
 
173
c      Time shift relative to data files <pfile0,sfile0> - format (hh.mm)
174
       read(10,*) timeshift
175
 
176
c      Read timeincrement between input files
177
       read(10,*) timeinc
178
 
179
c      Parameters for horizontal grid
180
       read(10,*) hmode
181
       if ( hmode.eq.'file' ) then                  ! from file
182
          read(10,*) hfile
183
       elseif ( hmode.eq.'line' ) then              ! along a line
184
          read(10,*) lon1,lon2,lat1,lat2,hn
185
       elseif ( hmode.eq.'box.eqd' ) then           ! box: 2d equidistant
186
          read(10,*) lon1,lon2,lat1,lat2,ds
187
       elseif ( hmode.eq.'box.grid' ) then          ! box: 2d grid
188
          read(10,*) lon1,lon2,lat1,lat2
189
       elseif ( hmode.eq.'point' ) then             ! single point
190
          read(10,*) lon1,lat1
191
       elseif ( hmode.eq.'shift' ) then             ! centre + shifted
192
          read(10,*) lon1,lat1,dlon,dlat
193
       elseif ( hmode.eq.'polygon.eqd' ) then       ! polygon: 2d equidistant
194
          read(10,*) hfile,ds
195
       elseif ( hmode.eq.'polygon.grid' ) then      ! polygon: 2d grid
196
          read(10,*) hfile
197
       elseif ( hmode.eq.'circle.eqd' ) then        ! circle: 2d equidistant
198
          read(10,*) lon1,lat1,radius,ds 
199
       elseif ( hmode.eq.'circle.grid' ) then       ! circle: 2d grid
200
          read(10,*) lon1,lat1,radius
201
       elseif ( hmode.eq.'region.eqd' ) then        ! region: 2d equidistant
202
          read(10,*) iregion,ds 
203
       elseif ( hmode.eq.'region.grid' ) then       ! iregion: 2d grid
204
          read(10,*) iregion
205
       elseif ( hmode.eq.'mask.grid' ) then         ! 0/1 masked - grid points
206
          read(10,*) hfile,maskname
207
       elseif ( hmode.eq.'mask.eqd' ) then         ! 0/1 masked - euqidistant points
208
          read(10,*) hfile,maskname,ds
209
       else
210
          print*,' ERROR: horizontal mode not supported ',trim(hmode)
211
          stop
212
       endif
213
 
214
c      Parameters for vertical grid
215
       read(10,*) vmode
216
       if ( vmode.eq.'file') then                   ! from file
217
          read(10,*) vfile
218
       elseif ( vmode.eq.'level' ) then             ! single level (explicit command)
219
          read(10,*) lev1                         
220
       elseif ( vmode.eq.'list') then               ! a list
221
          read(10,*) vn
222
          read(10,*) (levlist(i),i=1,vn)
223
       elseif ( vmode.eq.'profile') then            ! a profile
224
          read(10,*) lev1,lev2,vn
225
       elseif ( vmode.eq.'grid') then               ! grid points
226
          read(10,*) lev1,lev2
227
       else
228
          print*,' ERROR: vertical mode not supported ',trim(vmode)
229
          stop
230
       endif
231
 
232
c      Read units of vertical axis
233
       read(10,*) umode
234
       if ( ( umode.ne.'hPa'     ).and.
235
     >      ( umode.ne.'hPa,agl' ).and.   
236
     >      ( umode.ne.'K'       ).and.
237
     >      ( umode.ne.'PVU'     ).and.
238
     >      ( umode.ne.'INDEX'   )  ) 
239
     > then  
240
          print*,' ERROR: unit not supported ',trim(umode)
241
          stop
242
       endif
243
 
244
c     Read selection criterion (dummy read)
245
      read(10,*) selectstr
246
 
247
c     Read flag for 'no time check'
248
      read(10,*) timecheck 
249
 
250
c     Close parameter file
251
      close(10)
252
 
253
c     Decide which output format is used (1..4: trajectory format, -1: triple list)
254
      call mode_tra(oformat,ofile)
255
 
256
c     Decide whether all lat/lon/lev coordaintes are read from one file
257
      if ( (hmode.eq.'file').and.(vmode.eq.'nil') ) then
258
         hmode='file3'
259
      elseif ( (hmode.eq.'file').and.(vmode.ne.'nil') ) then
260
         hmode='file2'
261
      endif
262
 
263
c     Convert timeshift (hh.mm) into a fractional time shift
264
      call hhmm2frac(timeshift,tfrac)
265
      if (tfrac.gt.0.) then
266
         tfrac=tfrac/timeinc
267
      else
268
         tfrac=0.
269
      endif
270
 
271
c     If a mask file is provided, no time interpolation is allowed
272
      if ( (hmode.eq.'mask.grid').or.(hmode.eq.'mask.eqd') ) then
273
         if ( abs(tfrac).gt.eps ) then
274
            print*,' ERROR: no intermediate times allowed for ',
275
     >             trim(hmode)
276
            stop
277
         endif
278
      endif
279
 
280
c     Read the region coordinates if needed
281
      if ( (hmode.eq.'region.eqd' ).or.
282
     >     (hmode.eq.'region.grid') ) then
283
 
284
         open(10,file=regionf)
285
 
286
 50       read(10,*,end=51) string
287
 
288
          if ( string(1:1).ne.'#' ) then
289
             call regionsplit(string,i,xcorner,ycorner)
290
             if ( i.eq.iregion ) goto 52
291
          endif
292
 
293
          goto 50
294
 
295
 51      close(10)
296
 
297
         print*,' ERROR: region ',iregion,' not found on ',trim(regionf)
298
         stop
299
 
300
 52      continue
301
 
302
      endif
303
 
304
c     Write some status information
305
      print*,'---- INPUT PARAMETERS -----------------------------------'
306
      print*
307
      if ( timeshift.gt.0. ) then
308
         print*,'  P file                     : ',trim(pfile0),
309
     >                                            ' ',
310
     >                                            trim(pfile1)
311
         print*,'  S file                     : ',trim(sfile0),
312
     >                                            ' ',
313
     >                                            trim(sfile1)
314
      else
315
         print*,'  P file                     : ',trim(pfile0)
316
         print*,'  S file                     : ',trim(sfile0)
317
      endif
318
      print*,'  Output file                : ',trim(ofile) 
319
      print*
320
      if (oformat.eq.-1) then
321
         print*,'  Output format              : (lon,lat,lev)-list'
322
      else
323
         print*,'  Output format              : ',oformat
324
      endif
325
      print*
326
      print*,'  Reference time (year)      : ',reftime(1)
327
      print*,'                 (month)     : ',reftime(2)
328
      print*,'                 (day)       : ',reftime(3)
329
      print*,'                 (hour)      : ',reftime(4)
330
      print*,'                 (min)       : ',reftime(5)
331
      print*,'  Time range                 : ',reftime(6)
332
      print*
333
      print*,'  Time shift                 : ',timeshift,' + ',
334
     >                                         trim(pfile0)
335
      print*,'  Region file                : ',trim(regionf)
336
      print*
337
      print*,'  hmode                      : ',trim(hmode)
338
      if ( hmode.eq.'file2' ) then        
339
        print*,'      filename [lat/lon]      : ',trim(hfile)
340
      elseif ( hmode.eq.'file3' ) then        
341
        print*,'      filename [lat/lon/lev]  : ',trim(hfile)
342
      elseif ( hmode.eq.'line' ) then   
343
        write(*,'(a30,4f10.2,i4)') 
344
     >         '      lon1,lon2,lat1,lat2,n   : ',lon1,lon2,lat1,lat2,hn
345
      elseif ( hmode.eq.'box.eqd' ) then     
346
        write(*,'(a30,5f10.2)') 
347
     >         '      lon1,lon2,lat1,lat2,ds  : ',lon1,lon2,lat1,lat2,ds
348
      elseif ( hmode.eq.'box.grid' ) then    
349
        write(*,'(a30,4f10.2)') 
350
     >         '      lon1,lon2,lat1,lat2     : ',lon1,lon2,lat1,lat2 
351
      elseif ( hmode.eq.'point' ) then   
352
        print*,'      lon,lat                 : ',lon1,lat1
353
      elseif ( hmode.eq.'shift' ) then   
354
        write(*,'(a30,4f10.2)') 
355
     >         '      lon,lat,dlon,dlat       : ',lon1,lat1,dlon,dlat
356
      elseif ( hmode.eq.'polygon.eqd' ) then   
357
        write(*,'(a30,a10,f10.2)') 
358
     >         '      hfile, ds               : ',trim(hfile),ds
359
      elseif ( hmode.eq.'polygon.grid' ) then   
360
        write(*,'(a30,a10)') 
361
     >         '      hfile                   : ',trim(hfile)
362
      elseif ( hmode.eq.'circle.eqd' ) then   
363
        write(*,'(a30,4f10.2)') 
364
     >         '      lonc,latc,radius, ds    : ',lon1,lat1,radius,ds
365
      elseif ( hmode.eq.'circle.grid' ) then   
366
        write(*,'(a30,3f10.2)') 
367
     >         '      lonc,latc,radius        : ',lon1,lat1,radius
368
      elseif ( hmode.eq.'region.eqd' ) then   
369
        write(*,'(a30,i4,1f10.2)') 
370
     >         '      iregion, ds             : ',iregion,ds
371
        write(*,'(a30,4f10.2)')
372
     >         '      xcorner                 : ',(xcorner(i),i=1,4)
373
        write(*,'(a30,4f10.2)')
374
     >         '      ycorner                 : ',(ycorner(i),i=1,4)
375
      elseif ( hmode.eq.'region.grid' ) then   
376
        write(*,'(a30,i4)') 
377
     >         '      iregion                 : ',iregion
378
        write(*,'(a30,4f10.2)')
379
     >         '      xcorner                 : ',(xcorner(i),i=1,4)
380
        write(*,'(a30,4f10.2)')
381
     >         '      ycorner                 : ',(ycorner(i),i=1,4)
382
      elseif ( hmode.eq.'mask.eqd' ) then
383
        write(*,'(a30,a15,a15,f10.2)')
384
     >         '      hfile,variable,ds       : ',trim(hfile),
385
     >                                            trim(maskname),
386
     >                                            ds
387
      elseif ( hmode.eq.'mask.grid' ) then
388
        write(*,'(a30,a15,a15)')
389
     >         '      hfile,variable          : ',trim(hfile),
390
     >                                            trim(maskname)
391
      endif
392
 
393
      print*
394
      print*,'  vmode                      : ',trim(vmode)
395
      if ( vmode.eq.'file') then 
396
         print*,'      filename               : ',trim(vfile)
397
      elseif ( vmode.eq.'level' ) then 
398
         print*,'      level                  : ',lev1                         
399
      elseif ( vmode.eq.'list') then 
400
         print*,'      n                      : ',vn
401
         print*,'      level(i)               : ',(levlist(i),i=1,vn)
402
      elseif ( vmode.eq.'profile') then 
403
         print*,'      lev1,lev2,n            : ',lev1,lev2,vn
404
      elseif ( vmode.eq.'grid') then 
405
         print*,'      lev1,lev2              : ',lev1,lev2
406
      endif
407
 
408
      print* 
409
      print*,'  umode                      : ',trim(umode)
410
      print*
411
      print*,'  time check                 : ',trim(timecheck)
412
      print*
413
 
414
c     <mask.grid> and <mask.eqd>: split variable into name, operator,value
415
      if ( (hmode.eq.'mask.grid').or.(hmode.eq.'mask.eqd') ) then
416
          len = len_trim(maskname)
417
          do i=1,len_trim(maskname)-1
418
            ch1 = maskname(i:i)
419
            ch2 = maskname(i+1:i+1)
420
 
421
            if ( (ch1.eq.'<').and.(ch2.eq.'>') ) then
422
               read(maskname(i+2:len),*) maskvalue
423
               maskname = maskname(1:i-1)
424
               maskoper = 'ne'
425
               goto 90
426
            elseif ( (ch1.eq.'<').and.(ch2.eq.'=') ) then
427
               read(maskname(i+2:len),*) maskvalue
428
               maskname = maskname(1:i-1)
429
               maskoper = 'le'
430
               goto 90
431
            elseif ( (ch1.eq.'>').and.(ch2.eq.'=') ) then
432
               read(maskname(i+2:len),*) maskvalue
433
               maskname = maskname(1:i-1)
434
               maskoper = 'ge'
435
               goto 90
436
            elseif ( (ch1.eq.'=').and.(ch2.eq.'=') ) then
437
               read(maskname(i+2:len),*) maskvalue
438
               maskname = maskname(1:i-1)
439
               maskoper = 'eq'
440
               goto 90
441
            elseif ( ch1.eq.'=' ) then
442
               read(maskname(i+1:len),*) maskvalue
443
               maskname = maskname(1:i-1)
444
               maskoper = 'eq'
445
               goto 90
446
            elseif ( ch1.eq.'>' ) then
447
               read(maskname(i+1:len),*) maskvalue
448
               maskname = maskname(1:i-1)
449
               maskoper = 'gt'
450
               goto 90
451
            elseif ( ch1.eq.'<' ) then
452
               read(maskname(i+1:len),*) maskvalue
453
               maskname = maskname(1:i-1)
454
               maskoper = 'lt'
455
               goto 90
456
             endif
457
 
458
          enddo
459
      endif
460
 
461
  90  continue
462
 
463
c     ------------------------------------------------------------------
464
c     Read grid parameters from inital files
465
c     ------------------------------------------------------------------
466
 
467
c     Get the time of the first and second data file
468
      tstart = -timeshift                              ! Format hh.mm
469
      call hhmm2frac(tstart,frac)
470
      frac   = frac + timeinc          
471
      call frac2hhmm(frac,tend)                        ! Format hh.mm                  
472
 
473
c     Convert timeshift (hh.mm) into a fractional time shift
474
      tfrac=real(int(timeshift))+
475
     >          100.*(timeshift-real(int(timeshift)))/60.
476
      if (tfrac.gt.0.) then
477
         tfrac=tfrac/timeinc
478
      else
479
         tfrac=0.
480
      endif
481
 
482
c     Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,
483
c     pollon,pollat) The negative <-fid> of the file identifier is used 
484
c     as a flag for parameter retrieval
485
      varname  = 'U'
486
      nx       = 1
487
      ny       = 1
488
      nz       = 1
489
      call input_open (fid,pfile0)
490
      call input_grid (-fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
491
     >                 tstart,pollon,pollat,rd,rd,nz,rd,rd,rd,timecheck)
492
      call input_close(fid)
493
 
494
c     Check whether region coordinates are within the domain
495
      if ( (hmode.eq.'region.eqd' ).or.
496
     >     (hmode.eq.'region.grid') ) then
497
 
498
         do i=1,4
499
            if ( (xcorner(i).lt.xmin).or.
500
     >           (ycorner(i).lt.ymin).or.
501
     >           (xcorner(i).gt.xmax).or.
502
     >           (ycorner(i).gt.ymax) )
503
     >      then
504
               print*,' ERROR: region not included in data domain...'
505
               print*,'     ',trim(string)
506
               print*,'     ',(xcorner(j),j=1,4)
507
               print*,'     ',(ycorner(j),j=1,4)
508
               stop
509
            endif
510
 
511
         enddo
512
 
513
      endif
514
 
515
C     Check if the number of levels is too large
516
      if (nz.gt.nlevmax) goto 993
517
 
518
c     Allocate memory for 3d arrays: pressure, theta, pv
519
      allocate(pr(nx,ny,nz),stat=stat)
520
      if (stat.ne.0) print*,'*** error allocating array pr ***'
521
      allocate(prs(nx,ny),stat=stat)
522
      if (stat.ne.0) print*,'*** error allocating array prs **'
523
      allocate(th(nx,ny,nz),stat=stat)
524
      if (stat.ne.0) print*,'*** error allocating array th ***'
525
      allocate(ths(nx,ny),stat=stat)
526
      if (stat.ne.0) print*,'*** error allocating array ths **'
527
      allocate(pv(nx,ny,nz),stat=stat)
528
      if (stat.ne.0) print*,'*** error allocating array pv ***'
529
      allocate(pvs(nx,ny),stat=stat)
530
      if (stat.ne.0) print*,'*** error allocating array pvs **'
531
      allocate(in(nx,ny,nz),stat=stat)
532
      if (stat.ne.0) print*,'*** error allocating array in ***'
533
      allocate(mask(nx,ny),stat=stat)
534
      if (stat.ne.0) print*,'*** error allocating array mask ***'
535
 
536
c     Allocate memory for temporary arrays for time interpolation
537
      allocate(fld0(nx,ny,nz),stat=stat)
538
      if (stat.ne.0) print*,'*** error allocating array tmp0 ***'
539
      allocate(fld1(nx,ny,nz),stat=stat)
540
      if (stat.ne.0) print*,'*** error allocating array tmp1 ***'
541
      allocate(sfc0(nx,ny),stat=stat)
542
      if (stat.ne.0) print*,'*** error allocating array sfc0 ***'
543
      allocate(sfc1(nx,ny),stat=stat)
544
      if (stat.ne.0) print*,'*** error allocating array sfc1 ***'
545
 
546
c     ------ Index -----------------------------------------------------
547
 
548
c     Init the dummy array with vertical index
549
      do i=1,nx
550
         do j=1,ny
551
            do k=1,nz
552
               in(i,j,k) = real(k)
553
            enddo
554
         enddo
555
      enddo
556
 
557
c     ------ Pressure --------------------------------------------------
558
 
559
c     Read pressure from first data file (pfile0) on U-grid; we have to set
560
c     mdv explicitely, because it's not read from netCDF
561
      call input_open (fid,pfile0)
562
      varname='U'
563
      stagz  = -0.5
564
      call input_grid                                      
565
     >     (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
566
     >      tstart,pollon,pollat,fld0,sfc0,nz,ak,bk,stagz,timecheck)
567
      mdv = -999.99 
568
      call input_close(fid)
569
 
570
c     Read or set pressure for second data file (pfile1)
571
      if ( timeshift.ne.0.) then
572
         call input_open (fid,pfile1)
573
         varname='U'
574
         call input_grid                                      
575
     >        (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
576
     >        tend,pollon,pollat,fld1,sfc1,nz,ak,bk,stagz,timecheck)
577
         call input_close(fid)
578
      else
579
         do i=1,nx
580
            do j=1,ny
581
               do k=1,nz
582
                  fld1(i,j,k) = fld0(i,j,k)
583
               enddo
584
               sfc1(i,j) = sfc0(i,j)
585
            enddo
586
         enddo
587
      endif
588
 
589
c     Time interpolation to get the final pressure field
590
      do i=1,nx
591
         do j=1,ny
592
            do k=1,nz
593
               pr(i,j,k) = (1.-tfrac) * fld0(i,j,k) +
594
     >                     tfrac      * fld1(i,j,k)
595
            enddo
596
            prs(i,j) = (1.-tfrac) * sfc0(i,j) +
597
     >                 tfrac      * sfc1(i,j)
598
         enddo
599
      enddo
600
 
601
c     ------ Potential temperature -------------------------------------
602
 
603
       if ( (umode.eq.'K').or.(umode.eq.'PVU') ) then
604
 
605
c         Read potential temperature from first data file <sfile0>
606
          call input_open (fid,sfile0)
607
          varname='TH'                                      ! Theta                                      
608
          call input_wind
609
     >         (fid,varname,fld0,tstart,stagz,mdv,
610
     >         xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
611
          call input_close(fid)
612
 
613
c         Read or set potential temperature for second data file (sfile1)
614
          if ( timeshift.ne.0.) then
615
             call input_open (fid,sfile1)
616
             varname='TH' 
617
             call input_wind
618
     >            (fid,varname,fld1,tend,stagz,mdv,
619
     >            xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
620
             call input_close(fid)
621
          else
622
             do i=1,nx
623
                do j=1,ny
624
                   do k=1,nz
625
                      fld1(i,j,k) = fld0(i,j,k)
626
                   enddo
627
                enddo
628
             enddo
629
          endif
630
 
631
c         Time interpolation to get the final potential temperature field
632
          do i=1,nx
633
             do j=1,ny
634
                do k=1,nz
635
                   th(i,j,k) = (1.-tfrac) * fld0(i,j,k) +
636
     >                         tfrac      * fld1(i,j,k)
637
                enddo
638
             enddo
639
          enddo
640
 
641
c         Set the surface potential temperature
642
          do i=1,nx                            
643
             do j=1,ny
644
                ths(i,j)=th(i,j,1)
645
             enddo
646
          enddo
647
 
648
       endif
649
 
650
c     ------ Potential vorticity -----------------------------------------
651
 
652
       if ( (umode.eq.'PVU') ) then
653
 
654
c         Read potential vorticity from first data file <sfile0>
655
          call input_open (fid,sfile0)
656
          varname='PV'            
657
          call input_wind
658
     >         (fid,varname,fld0,tstart,stagz,mdv,
659
     >         xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
660
          call input_close(fid)
661
 
662
c         Read or set potential vorticity for second data file (sfile1)
663
          if ( timeshift.ne.0.) then
664
             call input_open (fid,sfile1)
665
             varname='PV' 
666
             call input_wind
667
     >            (fid,varname,fld1,tend,stagz,mdv,
668
     >            xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
669
             call input_close(fid)
670
          else
671
             do i=1,nx
672
                do j=1,ny
673
                   do k=1,nz
674
                      fld1(i,j,k) = fld0(i,j,k)
675
                   enddo
676
                enddo
677
             enddo
678
          endif
679
 
680
c         Time interpolation to get the final potential vorticity field
681
          do i=1,nx
682
             do j=1,ny
683
                do k=1,nz
684
                   pv(i,j,k) = (1.-tfrac) * fld0(i,j,k) +
685
     >                         tfrac      * fld1(i,j,k)
686
                enddo
687
             enddo
688
          enddo
689
 
690
c         Set the surface potential vorticity
691
          do i=1,nx                          
692
             do j=1,ny
693
                pvs(i,j)=pv(i,j,1)
694
             enddo
695
          enddo
696
       endif
697
 
698
c     ------ Load 0/1 label (mask) ------------------------------------
699
      if ( (hmode.eq.'mask.grid').or.(hmode.eq.'mask.eqd') ) then
700
 
701
c         Explicit netCDF calls - no consistency check for grids!
702
          ierr = NF90_OPEN(hfile,nf90_nowrite, cdfid)
703
          IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
704
          ierr = NF90_INQ_VARID(cdfid,maskname,varid)
705
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
706
          ierr = nf90_get_var(cdfid,varid,mask)
707
          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
708
          ierr = NF90_CLOSE(cdfid)
709
          IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
710
 
711
          nmask = 0
712
          do i=1,nx
713
            do j=1,ny
714
 
715
                if ( maskoper.eq.'eq' ) then
716
                   if ( abs(mask(i,j)-maskvalue).lt.eps ) then
717
                      nmask = nmask + 1
718
                      mask(i,j) = 1.
719
                   else
720
                      mask(i,j) = 0.
721
                   endif
722
                endif
723
 
724
                if ( maskoper.eq.'ne' ) then
725
                   if ( abs(mask(i,j)-maskvalue).gt.eps ) then
726
                      nmask = nmask + 1
727
                      mask(i,j) = 1.
728
                   else
729
                      mask(i,j) = 0.
730
                   endif
731
                endif
732
 
733
                if ( maskoper.eq.'gt' ) then
734
                   if ( mask(i,j).gt.maskvalue ) then
735
                      nmask = nmask + 1
736
                      mask(i,j) = 1.
737
                   else
738
                      mask(i,j) = 0.
739
                   endif
740
                endif
741
 
742
                if ( maskoper.eq.'lt' ) then
743
                   if ( mask(i,j).lt.maskvalue ) then
744
                      nmask = nmask + 1
745
                      mask(i,j) = 1.
746
                   else
747
                      mask(i,j) = 0.
748
                   endif
749
                endif
750
 
751
                if ( maskoper.eq.'ge' ) then
752
                   if ( mask(i,j).ge.maskvalue ) then
753
                      nmask = nmask + 1
754
                      mask(i,j) = 1.
755
                   else
756
                      mask(i,j) = 0.
757
                   endif
758
                endif
759
 
760
                if ( maskoper.eq.'le' ) then
761
                   if ( mask(i,j).le.maskvalue ) then
762
                      nmask = nmask + 1
763
                      mask(i,j) = 1.
764
                   else
765
                      mask(i,j) = 0.
766
                   endif
767
                endif
768
 
769
             enddo
770
          enddo
771
 
772
          print*
773
          print*,'  Mask read from file -> ',nmask,' 1-points'
774
          print*
775
 
776
      endif
777
 
778
c     Write some status information
779
      print*,'---- CONSTANT GRID PARAMETERS ---------------------------'
780
      print*
781
      print*,'  xmin,xmax        : ',xmin,xmax
782
      print*,'  ymin,ymax        : ',ymin,ymax
783
      print*,'  dx,dy            : ',dx,dy
784
      print*,'  pollon,pollat    : ',pollon,pollat
785
      print*,'  nx,ny,nz         : ',nx,ny,nz
786
      print*
787
      print*,'  Pressure loaded  : ',trim(pfile0),' ',trim(pfile1)
788
      if ( (umode.eq.'K').or.(umode.eq.'PVU') ) then
789
         print*,'  Theta loaded     : ',trim(sfile0),' ',trim(sfile1)
790
      endif
791
      if ( (umode.eq.'PVU') ) then
792
         print*,'  PV loaded        : ',trim(sfile0),' ',trim(sfile1)
793
      endif
794
      print*
795
 
796
c     ------------------------------------------------------------------
797
c     Determine the expanded list of starting coordinates
798
c     ------------------------------------------------------------------
799
 
800
c     Write some status information
801
      print*,'---- EXPAND LIST OF STARTING POSITIONS  -----------------'
802
      print*
803
 
804
c     ------ Read lat/lon/lev from <hfile> -----------------------------
805
      if ( hmode.eq.'file3' ) then
806
         start_n = 0
807
         open(10,file=hfile)
808
 100       continue
809
              start_n = start_n + 1
810
              read(10,*,end=101) start_lon(start_n),
811
     >                           start_lat(start_n),
812
     >                           start_lev(start_n)
813
              goto 100
814
 101       continue
815
           start_n = start_n - 1
816
         close(10)
817
         goto 400
818
      endif
819
 
820
c     ------ Get lat/lon (horizontal) coordinates ---------------------
821
 
822
c     Read lat/lon from <hfile> 
823
      if ( hmode.eq.'file2' ) then
824
         hn = 0
825
         open(10,file=hfile)
826
 200       continue
827
              hn = hn + 1
828
              read(10,*,end=201) lonlist(hn),
829
     >                           latlist(hn)
830
              goto 200
831
 201       continue
832
           hn = hn - 1
833
         close(10)
834
      endif
835
 
836
c     Get lat/lon along a line (linear in lat/lon space)
837
      if ( hmode.eq.'line' ) then
838
         do i=1,hn
839
            lonlist(i) = lon1 + real(i-1)/real(hn-1)*(lon2-lon1)
840
            latlist(i) = lat1 + real(i-1)/real(hn-1)*(lat2-lat1)
841
         enddo
842
      endif
843
 
844
c     Lat/lon box: equidistant
845
      if ( hmode.eq.'box.eqd' ) then
846
         hn  = 0
847
         lat = lat1
848
         do while ( lat.le.lat2 )      
849
           lon = lon1
850
           do while ( lon.le.lon2 ) 
851
             hn          = hn+1
852
             lonlist(hn) = lon
853
             latlist(hn) = lat
854
             lon         = lon+ds/(deltat*cos(pi180*lat))
855
           enddo
856
           lat = lat+ds/deltat
857
        enddo
858
      endif
859
 
860
c     Lat/lon box: grid
861
      if ( hmode.eq.'box.grid' ) then
862
         hn = 0
863
         do j=1,ny
864
            do i=1,nx
865
               lon = xmin + real(i-1) * dx
866
               lat = ymin + real(j-1) * dy
867
               if ( (lon.ge.lon1).and.(lon.le.lon2).and.
868
     >              (lat.ge.lat1).and.(lat.le.lat2) )
869
     >         then
870
                  hn          = hn+1
871
                  lonlist(hn) = lon
872
                  latlist(hn) = lat
873
               endif
874
            enddo
875
         enddo
876
      endif
877
 
878
c     Get single starting point
879
      if ( hmode.eq.'point' ) then
880
         hn          = 1
881
         lonlist(hn) = lon1
882
         latlist(hn) = lat1
883
      endif
884
 
885
c     Get shifted and central starting point
886
      if ( hmode.eq.'shift' ) then
887
         hn         = 5
888
         lonlist(1) = lon1
889
         latlist(1) = lat1
890
         lonlist(2) = lon1+dlon
891
         latlist(2) = lat1
892
         lonlist(3) = lon1-dlon
893
         latlist(3) = lat1
894
         lonlist(4) = lon1
895
         latlist(4) = lat1+dlat
896
         lonlist(5) = lon1
897
         latlist(5) = lat1-dlat
898
      endif
899
 
900
c     Lat/lon polygon: grid
901
      if ( hmode.eq.'polygon.grid' ) then
902
 
903
c        Read list of polygon coordinates
904
         pn = 0
905
         open(10,file=hfile)
906
           read(10,*) loninpoly,latinpoly
907
 210       continue
908
              pn = pn + 1
909
              read(10,*,end=211) lonpoly(pn),
910
     >                           latpoly(pn)
911
 
912
              print*,pn,lonpoly(pn),latpoly(pn)
913
 
914
              goto 210
915
 211       continue
916
           pn = pn - 1
917
         close(10)
918
 
919
c        Define the polygon boundaries
920
         call DefSPolyBndry(latpoly,lonpoly,pn,latinpoly,loninpoly)
921
 
922
c        Get the grid points inside the polygon
923
         hn = 0
924
         do j=1,ny
925
            do i=1,nx
926
               lon = xmin + real(i-1) * dx
927
               lat = ymin + real(j-1) * dy
928
 
929
               call LctPtRelBndry(lat,lon,flag)
930
 
931
               if ( (flag.eq.1).or.(flag.eq.2) ) then
932
                  hn          = hn+1
933
                  lonlist(hn) = lon
934
                  latlist(hn) = lat
935
               endif
936
 
937
            enddo
938
         enddo
939
 
940
      endif
941
 
942
c     Lat/lon polygon: equidistant
943
      if ( hmode.eq.'polygon.eqd' ) then
944
 
945
c        Read list of polygon coordinates
946
         pn = 0
947
 
948
         open(10,file=hfile)
949
           read(10,*) loninpoly,latinpoly
950
 220       continue
951
              pn = pn + 1
952
              read(10,*,end=221) lonpoly(pn),
953
     >                           latpoly(pn)
954
              goto 220
955
 221       continue
956
           pn = pn - 1
957
         close(10)
958
 
959
 
960
c        Define the polygon boundaries
961
         call DefSPolyBndry(latpoly,lonpoly,pn,latinpoly,loninpoly)
962
 
963
c        Get the grid points inside the polygon
964
         hn  = 0
965
         lat = -90.
966
         do while ( lat.le.90. )      
967
           lon = -180.
968
           do while ( lon.lt.180. )
969
 
970
              call LctPtRelBndry(lat,lon,flag)
971
 
972
               if ( (flag.eq.1).or.(flag.eq.2) ) then
973
                  hn          = hn+1
974
                  lonlist(hn) = lon
975
                  latlist(hn) = lat
976
 
977
               endif
978
 
979
               lon = lon+ds/(deltat*cos(pi180*lat))
980
           enddo
981
           lat = lat+ds/deltat
982
 
983
        enddo
984
 
985
      endif
986
 
987
c     Circle: equidistant
988
      if ( hmode.eq.'circle.eqd' ) then
989
         hn  = 0
990
         lat = ymin
991
         do while ( lat.le.ymax )      
992
           lon = xmin
993
           do while ( lon.le.xmax ) 
994
              dist = sdis(lon1,lat1,lon,lat)
995
              if ( dist.le.radius ) then
996
                 hn          = hn+1
997
                 lonlist(hn) = lon
998
                 latlist(hn) = lat
999
              endif
1000
              lon = lon+ds/(deltat*cos(pi180*lat))
1001
           enddo
1002
           lat = lat+ds/deltat
1003
        enddo
1004
      endif
1005
 
1006
c     Circle: grid
1007
      if ( hmode.eq.'circle.grid' ) then
1008
         hn = 0
1009
         do j=1,ny
1010
            do i=1,nx
1011
               lon = xmin + real(i-1) * dx
1012
               lat = ymin + real(j-1) * dy
1013
               dist = sdis(lon1,lat1,lon,lat)
1014
               if ( dist.le.radius ) then
1015
                  hn          = hn+1
1016
                  lonlist(hn) = lon
1017
                  latlist(hn) = lat
1018
               endif
1019
            enddo
1020
         enddo
1021
 
1022
      endif
1023
 
1024
c     Region: equidistant
1025
      if ( hmode.eq.'region.eqd' ) then
1026
         hn  = 0
1027
         lat = ymin
1028
         do while ( lat.le.ymax )      
1029
           lon = xmin
1030
           do while ( lon.le.xmax ) 
1031
              flag = inregion(lon,lat,xcorner,ycorner)
1032
              if ( flag.eq.1 ) then
1033
                 hn          = hn+1
1034
                 lonlist(hn) = lon
1035
                 latlist(hn) = lat
1036
              endif
1037
              lon = lon+ds/(deltat*cos(pi180*lat))
1038
           enddo
1039
           lat = lat+ds/deltat
1040
        enddo
1041
      endif
1042
 
1043
c     Region: grid
1044
      if ( hmode.eq.'region.grid' ) then
1045
         hn = 0
1046
         do j=1,ny
1047
            do i=1,nx
1048
               lon  = xmin + real(i-1) * dx
1049
               lat  = ymin + real(j-1) * dy
1050
               flag = inregion(lon,lat,xcorner,ycorner)
1051
               if ( flag.eq.1 ) then
1052
                  hn          = hn+1
1053
                  lonlist(hn) = lon
1054
                  latlist(hn) = lat
1055
               endif
1056
            enddo
1057
         enddo
1058
      endif
1059
 
1060
c     mask: grid
1061
      if ( hmode.eq.'mask.grid' ) then
1062
         hn  = 0
1063
         do i=1,nx
1064
          do j=1,ny
1065
            if ( mask(i,j).gt.0.5 ) then
1066
               hn = hn +1
1067
               lonlist(hn) = xmin + real(i-1) * dx
1068
               latlist(hn) = ymin + real(j-1) * dy
1069
            endif
1070
          enddo
1071
         enddo
1072
      endif
1073
 
1074
c     mask: equidistant
1075
      if ( hmode.eq.'mask.eqd' ) then
1076
         hn  = 0
1077
         lat = -90.+dy
1078
         do while ( lat.le.ymax )
1079
           lon = -180.
1080
           do while ( lon.le.xmax )
1081
              indx = nint( (lon-xmin)/dx + 1. )
1082
              indy = nint( (lat-ymin)/dy + 1. )
1083
              if ( indx.lt. 1 ) indx =  1
1084
              if ( indy.lt. 1 ) indy =  1
1085
              if ( indx.gt.nx ) indx = nx
1086
              if ( indy.gt.ny ) indy = ny
1087
              if ( mask(indx,indy).gt.0.5 ) then
1088
                 hn          = hn+1
1089
                 lonlist(hn) = lon
1090
                 latlist(hn) = lat
1091
              endif
1092
              lon = lon+ds/(deltat*cos(pi180*lat))
1093
           enddo
1094
           lat = lat+ds/deltat
1095
        enddo
1096
      endif
1097
 
1098
 
1099
c     ------ Get lev (vertical) coordinates -------------------------
1100
 
1101
c     Read level list from file
1102
      if ( vmode.eq.'file' ) then
1103
         vn = 0
1104
         open(10,file=vfile)
1105
 300       continue
1106
              vn = vn + 1
1107
              read(10,*,end=301) levlist(vn)
1108
              goto 300
1109
 301       continue
1110
           vn = vn - 1
1111
         close(10)
1112
      endif
1113
 
1114
c     Get single starting level
1115
      if ( vmode.eq.'level' ) then
1116
         vn          = 1
1117
         levlist(vn) = lev1
1118
      endif
1119
 
1120
c     Get level profile
1121
      if ( vmode.eq.'profile' ) then
1122
         do i=1,vn
1123
            levlist(i) = lev1 + real(i-1)/real(vn-1)*(lev2-lev1)
1124
         enddo
1125
      endif
1126
 
1127
c     Get all grid points in a layer: at the moment set the list of levels to 
1128
c     all indices from 1 to nz; later the correct subset of indices will be chosen
1129
      if ( vmode.eq.'grid' ) then
1130
         vn = nz
1131
         do i=1,vn
1132
            levlist(i) = real(i)
1133
         enddo
1134
         umode_save = umode
1135
         umode      = 'INDEX'
1136
 
1137
      endif
1138
 
1139
c     ------ Compile the complete list of starting positions  ------
1140
 
1141
c     Get all starting points in specified vertical coordinate system
1142
      start_n = 0
1143
      do i=1,vn
1144
         do j=1,hn
1145
 
1146
            start_n = start_n + 1
1147
            if ( start_n.ge.nmax ) then
1148
               print*,' ERROR: maximum number of start points exceeded'
1149
               stop
1150
            endif
1151
            start_lon(start_n) = lonlist(j)
1152
            start_lat(start_n) = latlist(j)
1153
            start_lev(start_n) = levlist(i)
1154
 
1155
         enddo
1156
      enddo
1157
 
1158
c     ------ Exit point of this section
1159
 400  continue
1160
 
1161
c     Write status information
1162
      print*,'  # expanded points : ', start_n
1163
      print*
1164
 
1165
c     ------------------------------------------------------------------
1166
c     Transform starting levels into pressure
1167
c     ------------------------------------------------------------------
1168
 
1169
c     Write some status information
1170
      print*,'---- STARTING POSITIONS ---------------------------------'
1171
      print*
1172
 
1173
c     Vertical mode <hPa,asl> or simply <hPa>
1174
      if ( (umode.eq.'hPa,asl').or.(umode.eq.'hPa') ) then
1175
 
1176
         do i=1,start_n
1177
            start_pre(i) = start_lev(i)
1178
         enddo
1179
 
1180
c     Vertical mode <hPa,agl>
1181
      elseif ( umode.eq.'hPa,agl' ) then
1182
 
1183
         do i=1,start_n
1184
            call get_index3(rid,rjd,rkd,start_lon(i),start_lat(i),1050.,
1185
     >                      3,pr,prs,nx,ny,nz,xmin,ymin,dx,dy)
1186
            tmp1 = int_index3 (prs,nx,ny,1,rid,rjd,1,mdv)
1187
            start_pre(i) = tmp1 - start_lev(i) 
1188
         enddo
1189
 
1190
c     Vertical mode <K>
1191
      elseif ( umode.eq.'K' ) then
1192
 
1193
         do i=1,start_n
1194
            call get_index3(rid,rjd,rkd,start_lon(i),start_lat(i),
1195
     >                start_lev(i),1,th,ths,nx,ny,nz,xmin,ymin,dx,dy)
1196
            tmp1 = int_index3 (pr,nx,ny,nz,rid,rjd,rkd,mdv)
1197
            start_pre(i) = tmp1
1198
         enddo
1199
 
1200
c     Vertical mode <PVU>
1201
      elseif ( umode.eq.'PVU' ) then
1202
 
1203
         do i=1,start_n
1204
            call get_index3(rid,rjd,rkd,start_lon(i),start_lat(i),
1205
     >               start_lev(i),2,pv,pvs,nx,ny,nz,xmin,ymin,dx,dy)
1206
            tmp1 = int_index3 (pr,nx,ny,nz,rid,rjd,rkd,mdv)
1207
            start_pre(i) = tmp1
1208
         enddo
1209
 
1210
c     Vertical mode <INDEX>
1211
      elseif ( umode.eq.'INDEX' ) then
1212
 
1213
         do i=1,start_n
1214
            call get_index3(rid,rjd,rkd,start_lon(i),start_lat(i),
1215
     >               1050.,2,pv,pvs,nx,ny,nz,xmin,ymin,dx,dy)
1216
            rkd = start_lev(i)
1217
            tmp1 = int_index3 (pr,nx,ny,nz,rid,rjd,rkd,mdv)
1218
            start_pre(i) = tmp1
1219
         enddo
1220
 
1221
      endif
1222
 
1223
c     ------------------------------------------------------------------
1224
c     Remove invalid points from the list
1225
c     ------------------------------------------------------------------
1226
 
1227
c     Select the correct subset if <vmode=grid>: starting points outside the layer
1228
c     will receive a <mdv> vertical pressure and will be removed
1229
 
1230
      if ( vmode.eq.'grid' ) then
1231
 
1232
         do i=1,start_n
1233
 
1234
c           Get the pressure at the grid point
1235
            if ( ( umode_save.eq.'hPa'      ).or.
1236
     >            (umode_save.eq.'hPa,asl') ) 
1237
     >      then
1238
 
1239
               call get_index3(rid,rjd,rkd,start_lon(i),start_lat(i),
1240
     >                         start_pre(i),3,pr,prs,nx,ny,nz,xmin,
1241
     >                         ymin,dx,dy)
1242
               tmp1 = int_index3 (pr,nx,ny,nz,rid,rjd,rkd,mdv)
1243
 
1244
c           Get pressure AGL at grid point
1245
            elseif ( umode_save.eq.'hPa,agl' ) then
1246
 
1247
               call get_index3(rid,rjd,rkd,start_lon(i),
1248
     >                         start_lat(i),start_pre(i),3,pr,prs,
1249
     >                         nx,ny,nz,xmin,ymin,dx,dy)
1250
               tmp1 = int_index3 (pr,nx,ny,nz,rid,rjd,rkd,mdv)
1251
               call get_index3(rid,rjd,rkd,start_lon(i),
1252
     >                         start_lat(i),1050.,3,pr,prs,nx,ny,
1253
     >                         nz,xmin,ymin,dx,dy)
1254
               tmp2 = int_index3 (prs,nx,ny,1,rid,rjd,1,mdv)
1255
               tmp1 = tmp2 - tmp1
1256
 
1257
c           Get potential temperature at grid point
1258
            elseif ( umode_save.eq.'K' ) then
1259
 
1260
               call get_index3(rid,rjd,rkd,start_lon(i),start_lat(i),
1261
     >                         start_pre(i),3,pr,prs,nx,ny,nz,
1262
     >                         xmin,ymin,dx,dy)
1263
               tmp1 = int_index3 (th,nx,ny,nz,rid,rjd,rkd,mdv)
1264
 
1265
c           Get potential vorticity at the grid point
1266
            elseif ( umode_save.eq.'PVU' ) then
1267
 
1268
               call get_index3(rid,rjd,rkd,start_lon(i),start_lat(i),
1269
     >                         start_pre(i),3,pr,prs,nx,ny,nz,xmin,
1270
     >                         ymin,dx,dy)
1271
               tmp1 = int_index3 (pv,nx,ny,nz,rid,rjd,rkd,mdv)
1272
 
1273
c           Get vertical index at the grid point
1274
            elseif ( umode_save.eq.'INDEX' ) then
1275
 
1276
              call get_index3(rid,rjd,rkd,start_lon(i),start_lat(i),
1277
     >                        start_pre(i),3,pr,prs,nx,ny,nz,
1278
     >                        xmin,ymin,dx,dy)
1279
              tmp1 = int_index3 (in,nx,ny,nz,rid,rjd,rkd,mdv)
1280
 
1281
           endif
1282
 
1283
c          Remove points outside layer
1284
           if ( ( tmp1.lt.lev1).or.(tmp1.gt.lev2) ) then
1285
              start_pre(i) = mdv
1286
           endif
1287
 
1288
         enddo         
1289
 
1290
      endif
1291
 
1292
c     Check whether the starting levels are valid (in data domain)
1293
      do i=1,start_n
1294
 
1295
         call get_index3(rid,rjd,rkd,start_lon(i),start_lat(i),1050.,
1296
     >                   3,pr,prs,nx,ny,nz,xmin,ymin,dx,dy)
1297
c         print*,'prs',prs
1298
         tmp1 = int_index3 (prs,nx,ny, 1,rid,rjd,real( 1),mdv)           ! Surface
1299
         tmp2 = int_index3 (pr ,nx,ny,nz,rid,rjd,real(nz),mdv)           ! Top of domain
1300
 
1301
         if ( (start_pre(i).gt.tmp1).or.
1302
     >        (start_pre(i).lt.tmp2).or.
1303
     >        (start_lon(i).lt.xmin).or. 
1304
     >        (start_lon(i).gt.xmax).or. 
1305
     >        (start_lat(i).lt.ymin).or.
1306
     >        (start_lat(i).gt.ymax) )
1307
     >   then
1308
            start_pre(i) = mdv
1309
         endif
1310
 
1311
      enddo
1312
 
1313
c     Mark all starting points outside the domain
1314
      n_outside = 0
1315
      do i=1,start_n
1316
 
1317
         start_valid(i) = 1
1318
         if ( abs(start_pre(i)-mdv).lt.eps ) then
1319
 
1320
            start_valid(i) = 0
1321
            if ( vmode.ne.'grid' ) then
1322
               if ( n_outside.lt.10 ) then
1323
                  print*,'  Outside ', 
1324
     >                 start_lon(i),start_lat(i),start_lev(i)
1325
                  n_outside = n_outside + 1
1326
               elseif ( n_outside.eq.10 ) then
1327
                  print*,'  Outside: more than 10 starting points'
1328
                  n_outside = n_outside + 1
1329
               else
1330
                  n_outside = n_outside + 1
1331
               endif
1332
            endif
1333
 
1334
         endif    
1335
 
1336
      enddo 
1337
 
1338
 
1339
c     Only keep valid starting points 
1340
      j = 0
1341
      do i=1,start_n
1342
         if ( start_valid(i).eq.1 ) then
1343
             j = j + 1
1344
             start_lon(j) = start_lon(i)
1345
             start_lat(j) = start_lat(i)
1346
             start_pre(j) = start_pre(i)
1347
             start_lev(j) = start_lev(i)
1348
          endif
1349
      enddo
1350
      start_n = j
1351
 
1352
c     Write some status information
1353
      latmin = start_lat(1)
1354
      latmax = start_lat(1)
1355
      lonmin = start_lon(1)
1356
      lonmax = start_lon(1)
1357
      premin = start_pre(1)
1358
      premax = start_pre(1)
1359
      do i=1,start_n
1360
         if (start_lat(i).lt.latmin) latmin = start_lat(i)
1361
         if (start_lat(i).gt.latmax) latmax = start_lat(i)
1362
         if (start_lon(i).lt.lonmin) lonmin = start_lon(i)
1363
         if (start_lon(i).gt.lonmax) lonmax = start_lon(i)
1364
         if (start_pre(i).lt.premin) premin = start_pre(i)
1365
         if (start_pre(i).gt.premax) premax = start_pre(i)
1366
      enddo
1367
      print*,'  min(lat),max(lat) : ', latmin,latmax
1368
      print*,'  min(lon),max(lon) : ', lonmin,lonmax
1369
      print*,'  min(pre),max(pre) : ', premin,premax
1370
      print*
1371
      print*,'  # starting points : ', start_n
1372
      print*,'  # Ouside domain   : ', n_outside
1373
      print*
1374
 
1375
 
1376
c     ------------------------------------------------------------------
1377
c     Write starting positions to output file
1378
c     ------------------------------------------------------------------
1379
 
1380
c     Output as a trajectory file (with only one time == 0)
1381
      if (oformat.ne.-1) then
1382
 
1383
         allocate(tra(start_n,1,5),stat=stat)
1384
 
1385
         vars(1)  ='time'
1386
         vars(2)  ='lon'
1387
         vars(3)  ='lat'
1388
         vars(4)  ='p'
1389
         vars(5)  ='level'
1390
 
1391
         call wopen_tra(fid,ofile,start_n,1,5,reftime,vars,oformat)
1392
 
1393
         do i=1,start_n
1394
            tra(i,1,1) = 0.
1395
            tra(i,1,2) = start_lon(i)
1396
            tra(i,1,3) = start_lat(i)
1397
            tra(i,1,4) = start_pre(i)
1398
            tra(i,1,5) = start_lev(i)
1399
         enddo
1400
         call write_tra(fid,tra,start_n,1,5,oformat)
1401
 
1402
         call close_tra(fid,oformat)
1403
 
1404
c     Output as a triple list (corresponding to <startf> file)
1405
      else
1406
 
1407
         fid = 10
1408
         open(fid,file=ofile)
1409
          do i=1,start_n
1410
             write(fid,'(3f10.3)') start_lon(i),start_lat(i),
1411
     >                             start_pre(i) 
1412
          enddo
1413
         close(fid)
1414
 
1415
      endif
1416
 
1417
c     Write some status information, and end of program message
1418
      print*  
1419
      print*,'---- STATUS INFORMATION --------------------------------'
1420
      print*
1421
      print*,'ok'
1422
      print*
1423
      print*,'       *** END OF PROGRAM CREATE_STARTF ***'
1424
      print*,'========================================================='
1425
 
1426
c     ------------------------------------------------------------------
1427
c     Exception handling
1428
c     ------------------------------------------------------------------
1429
 
1430
      stop
1431
 
1432
 993  write(*,*) '*** ERROR: problems with array size'
1433
      call exit(1)
1434
 
1435
      end
1436
 
1437
c     --------------------------------------------------------------------------
1438
c     Split a region string and get corners of the domain
1439
c     --------------------------------------------------------------------------
1440
 
1441
      subroutine regionsplit(string,iregion,xcorner,ycorner)
1442
 
1443
c     The region string comes either as <lonw,lone,lats,latn> or as <lon1,lat1,
1444
c     lon2,lat2,lon3,lat3,lon4,lat4>: split it into ints components and get the
1445
c     four coordinates for the region
1446
 
1447
      implicit none
1448
 
1449
c     Declaration of subroutine parameters
1450
      character*80    string
1451
      real            xcorner(4),ycorner(4)
1452
      integer         iregion
1453
 
1454
c     Local variables
1455
      integer         i,n
1456
      integer         il,ir
1457
      real            subfloat (80)
1458
      integer         stat
1459
      integer         len
1460
 
1461
c     ------- Split the string
1462
      i    = 1
1463
      n    = 0
1464
      stat = 0
1465
      il   = 1
1466
      len  = len_trim(string)
1467
 
1468
 100  continue
1469
 
1470
c     Find start of a substring
1471
      do while ( stat.eq.0 )
1472
         if ( string(i:i).ne.' ' ) then
1473
            stat = 1
1474
            il   = i
1475
         else
1476
            i = i + 1
1477
         endif
1478
      enddo
1479
 
1480
c     Find end of substring
1481
      do while ( stat.eq.1 )         
1482
         if ( ( string(i:i).eq.' ' ) .or. ( i.eq.len ) ) then
1483
            stat = 2
1484
            ir   = i
1485
         else
1486
            i    = i + 1
1487
         endif
1488
      enddo
1489
 
1490
c     Convert the substring into a number
1491
      if ( stat.eq.2 ) then
1492
         n = n + 1
1493
         read(string(il:ir),*) subfloat(n)
1494
         stat = 0
1495
      endif
1496
 
1497
      if ( i.lt.len ) goto 100
1498
 
1499
 
1500
c     -------- Get the region number
1501
 
1502
      iregion = nint(subfloat(1))
1503
 
1504
c     -------- Get the corners of the region
1505
 
1506
      if ( n.eq.5 ) then     ! lonw(2),lone(3),lats(4),latn(5)
1507
 
1508
         xcorner(1) = subfloat(2)
1509
         ycorner(1) = subfloat(4)
1510
 
1511
         xcorner(2) = subfloat(3)
1512
         ycorner(2) = subfloat(4)
1513
 
1514
         xcorner(3) = subfloat(3)
1515
         ycorner(3) = subfloat(5)
1516
 
1517
         xcorner(4) = subfloat(2)
1518
         ycorner(4) = subfloat(5)
1519
 
1520
      elseif ( n.eq.9 ) then     ! lon1,lat1,lon2,lat2,lon3,lon4,lat4
1521
 
1522
         xcorner(1) = subfloat(2)
1523
         ycorner(1) = subfloat(3)
1524
 
1525
         xcorner(2) = subfloat(4)
1526
         ycorner(2) = subfloat(5)
1527
 
1528
         xcorner(3) = subfloat(6)
1529
         ycorner(3) = subfloat(7)
1530
 
1531
         xcorner(4) = subfloat(8)
1532
         ycorner(4) = subfloat(9)
1533
 
1534
      else
1535
 
1536
         print*,' ERROR: invalid region specification '
1537
         print*,'     ',trim(string)
1538
         stop
1539
 
1540
      endif
1541
 
1542
 
1543
      end
1544
 
1545
c     --------------------------------------------------------------------------
1546
c     Decide whether lat/lon point is in or out of region
1547
c     --------------------------------------------------------------------------
1548
 
1549
      integer function inregion (lon,lat,xcorner,ycorner)
1550
 
1551
c     Decide whether point (lon/lat) is in the region specified by <xcorner(1..4),
1552
c     ycorner(1..4).
1553
 
1554
      implicit none
1555
 
1556
c     Declaration of subroutine parameters
1557
      real    lon,lat
1558
      real    xcorner(4),ycorner(4)
1559
 
1560
c     Local variables
1561
      integer flag
1562
      real    xmin,xmax,ymin,ymax
1563
      integer i
1564
 
1565
c     Reset the flag
1566
      flag = 0
1567
 
1568
c     Set some boundaries
1569
      xmax = xcorner(1)
1570
      xmin = xcorner(1)
1571
      ymax = ycorner(1)
1572
      ymin = ycorner(1)
1573
      do i=2,4
1574
        if (xcorner(i).lt.xmin) xmin = xcorner(i)
1575
        if (xcorner(i).gt.xmax) xmax = xcorner(i)
1576
        if (ycorner(i).lt.ymin) ymin = ycorner(i)
1577
        if (ycorner(i).gt.ymax) ymax = ycorner(i)
1578
      enddo
1579
 
1580
c     Do the tests - set flag=1 if all tests pased
1581
      if (lon.lt.xmin) goto 970
1582
      if (lon.gt.xmax) goto 970
1583
      if (lat.lt.ymin) goto 970
1584
      if (lat.gt.ymax) goto 970
1585
 
1586
      if ((lon-xcorner(1))*(ycorner(2)-ycorner(1))-
1587
     >    (lat-ycorner(1))*(xcorner(2)-xcorner(1)).gt.0.) goto 970
1588
      if ((lon-xcorner(2))*(ycorner(3)-ycorner(2))-
1589
     >    (lat-ycorner(2))*(xcorner(3)-xcorner(2)).gt.0.) goto 970
1590
      if ((lon-xcorner(3))*(ycorner(4)-ycorner(3))-
1591
     >    (lat-ycorner(3))*(xcorner(4)-xcorner(3)).gt.0.) goto 970
1592
      if ((lon-xcorner(4))*(ycorner(1)-ycorner(4))-
1593
     >    (lat-ycorner(4))*(xcorner(1)-xcorner(4)).gt.0.) goto 970
1594
 
1595
      flag = 1
1596
 
1597
c     Return the value
1598
 970  continue
1599
 
1600
      inregion = flag
1601
 
1602
      return
1603
 
1604
      end
1605
 
1606
c     --------------------------------------------------------------------------
1607
c     Spherical distance between lat/lon points                                                       
1608
c     --------------------------------------------------------------------------
1609
 
1610
      real function sdis(xp,yp,xq,yq)
1611
c
1612
c     calculates spherical distance (in km) between two points given
1613
c     by their spherical coordinates (xp,yp) and (xq,yq), respectively.
1614
c
1615
      real      re
1616
      parameter (re=6370.)
1617
      real      pi180
1618
      parameter (pi180=3.14159/180.)
1619
      real      xp,yp,xq,yq,arg
1620
 
1621
      arg=sin(pi180*yp)*sin(pi180*yq)+
1622
     >    cos(pi180*yp)*cos(pi180*yq)*cos(pi180*(xp-xq))
1623
      if (arg.lt.-1.) arg=-1.
1624
      if (arg.gt.1.) arg=1.
1625
 
1626
      sdis=re*acos(arg)
1627
 
1628
      end
1629
 
1630
 
1631
c     ****************************************************************
1632
c     * Given some spherical polygon S and some point X known to be  *
1633
c     * located inside S, these routines will determine if an arbit- *
1634
c     * -rary point P lies inside S, outside S, or on its boundary.  *
1635
c     * The calling program must first call DefSPolyBndry to define  *
1636
c     * the boundary of S and the point X. Any subsequent call to    *
1637
c     * subroutine LctPtRelBndry will determine if some point P lies *
1638
c     * inside or outside S, or on its boundary. (Usually            *
1639
c     * DefSPolyBndry is called once, then LctPrRelBndry is called   *
1640
c     * many times).                                                 *
1641
c     *                                                              * 
1642
c     * REFERENCE:            Bevis, M. and Chatelain, J.-L. (1989)  * 
1643
c     *                       Maflaematical Geology, vol 21.         *
1644
c     * VERSION 1.0                                                  *
1645
c     ****************************************************************
1646
 
1647
      Subroutine DefSPolyBndry(vlat,vlon,nv,xlat, xlon)
1648
 
1649
c     ****************************************************************
1650
c     * This mmn entry point is used m define ~e spheric~ polygon S  *
1651
c     * and the point X.                                             *
1652
c     * ARGUMENTS:                                                   *
1653
c     * vlat,vlon (sent) ... vectors containing the latitude and     * 
1654
c     *                      longitude of each vertex of the         *
1655
c     *                      spherical polygon S. The ith.vertex is  *
1656
c     *                      located at [vlat(i),vlon(i)].           *
1657
c     * nv        (sent) ... the number of vertices and sides in the *
1658
c     *                      spherical polygon S                     *
1659
c     * xlat,xlon (sent) ... latitude and longitude of some point X  *
1660
c     *                      located inside S. X must not be located *
1661
c     *                      on any great circle that includes two   *
1662
c     *                      vertices of S.                          *
1663
c     *                                                              *
1664
c     * UNITS AND SIGN CONVENTION:                                   *
1665
c     *  Latitudes and longitudes are specified in degrees.          *
1666
c     *  Latitudes are positive to the north and negative to the     *
1667
c     *  south.                                                      *
1668
c     *  Longitudes are positive to the east and negative to the     *
1669
c     *  west.                                                       *
1670
c     *                                                              * 
1671
c     * VERTEX ENUMERATION:                                          * 
1672
c     * The vertices of S should be numbered sequentially around the *
1673
c     * border of the spherical polygon. Vertex 1 lies between vertex*
1674
c     * nv and vertex 2. Neighbouring vertices must be seperated by  *
1675
c     * less than 180 degrees. (In order to generate a polygon side  *
1676
c     * whose arc length equals or exceeds 180 degrees simply        *
1677
c     * introduce an additional (pseudo)vertex). Having chosen       *
1678
c     * vertex 1, the user may number the remaining vertices in      *
1679
c     * either direction. However if the user wishes to use the      *
1680
c     * subroutine SPA to determine the area of the polygon S (Bevis *
1681
c     * & Cambareri, 1987, Math. Geol., v.19, p. 335-346) then he or *
1682
c     * she must follow the convention whereby in moving around the  *
1683
c     * polygon border in the direction of increasing vertex number  *
1684
c     * clockwise bends occur at salient vertices. A vertex is       *
1685
c     * salient if the interior angle is less than 180 degrees.      *
1686
c     * (In the case of a convex polygon this convention implies     *
1687
c     * that vertices are numbered in clockwise sequence).           *
1688
c     ****************************************************************
1689
 
1690
      implicit none
1691
 
1692
      integer mxnv,nv
1693
 
1694
c     ----------------------------------------------------------------
1695
c     Edit next statement to increase maximum number of vertices that 
1696
c     may be used to define the spherical polygon S               
1697
c     The value of parameter mxnv in subroutine LctPtRelBndry must match
1698
c     that of parameter mxnv in this subroutine, as assigned above.
1699
c     ----------------------------------------------------------------
1700
      parameter (mxnv=500)
1701
 
1702
      real  vlat(nv),vlon(nv),xlat,xlon,dellon
1703
      real  tlonv(mxnv),vlat_c(mxnv),vlon_c(mxnv),xlat_c,xlon_c
1704
      integer i,ibndry,nv_c,ip
1705
 
1706
      data ibndry/0/
1707
 
1708
      common/spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry
1709
 
1710
      if (nv.gt.mxnv) then
1711
         print *,'nv exceeds maximum allowed value'
1712
         print *,'adjust parameter mxnv in subroutine DefSPolyBndry'
1713
         stop
1714
      endif
1715
 
1716
      ibndry=1                  ! boundary defined at least once (flag)
1717
      nv_c=nv                   ! copy for named common
1718
      xlat_c=xlat               ! . . . .
1719
      xlon_c=xlon               !
1720
 
1721
      do i=1,nv
1722
         vlat_c(i)=vlat(i)      ! "
1723
         vlon_c(i)=vlon(i)      !
1724
 
1725
         call TrnsfmLon(xlat,xlon,vlat(i),vlon(i),tlonv(i))
1726
 
1727
         if (i.gt.1) then
1728
            ip=i-1
1729
         else
1730
            ip=nv
1731
         endif
1732
 
1733
         if ((vlat(i).eq.vlat(ip)).and.(vlon(i).eq.vlon(ip))) then
1734
            print *,'DefSPolyBndry detects user error:'
1735
            print *,'vertices ',i,' and ',ip,' are not distinct'
1736
            print*,'lat ',i,ip,vlat(i),vlat(ip)
1737
            print*,'lon ',i,ip,vlon(i),vlon(ip)            
1738
            stop
1739
         endif
1740
 
1741
         if (tlonv(i).eq.tlonv(ip)) then
1742
            print *,'DefSPolyBndry detects user error:'
1743
            print *,'vertices ',i,' & ',ip,' on same gt. circle as X'
1744
            stop
1745
         endif
1746
 
1747
         if (vlat(i).eq.(-vlat(ip))) then
1748
            dellon=vlon(i)-vlon(ip)
1749
            if (dellon.gt.+180.) dellon=dellon-360.
1750
            if (dellon.lt.-180.) dellon=dellon-360.
1751
            if ((dellon.eq.+180.0).or.(dellon.eq.-180.0)) then
1752
               print *,'DefSPolyBndry detects user error:'
1753
               print *,'vertices ',i,' and ',ip,' are antipodal'
1754
               stop
1755
            endif
1756
         endif
1757
      enddo
1758
 
1759
      return
1760
 
1761
      end
1762
 
1763
 
1764
c     ****************************************************************
1765
 
1766
      Subroutine LctPtRelBndry(plat,plon,location)
1767
 
1768
c     ****************************************************************
1769
 
1770
c     ****************************************************************
1771
c     * This routine is used to see if some point P is located       *
1772
c     * inside, outside or on the boundary of the spherical polygon  *
1773
c     * S previously defined by a call to subroutine DefSPolyBndry.  *
1774
c     * There is a single restriction on point P: it must not be     *
1775
c     * antipodal to the point X defined in the call to DefSPolyBndry*
1776
c     * (ie.P and X cannot be seperated by exactly 180 degrees).     *
1777
c     * ARGUMENTS:                                                   *  
1778
c     * plat,plon (sent)... the latitude and longitude of point P    *
1779
c     * location (returned)... specifies the location of P:          *
1780
c     *                        location=0 implies P is outside of S  *
1781
c     *                        location=1 implies P is inside of S   *
1782
c     *                        location=2 implies P on boundary of S *
1783
c     *                        location=3 implies user error (P is   *
1784
c     *                                     antipodal to X)          *
1785
c     * UNFfS AND SIGN CONVENTION:                                   * 
1786
c     *  Latitudes and longitudes are specified in degrees.          *
1787
c     *  Latitudes are positive to the north and negative to the     *
1788
c     *  south.                                                      *    
1789
c     *  Longitudes are positive to the east and negative to the     *
1790
c     *  west.                                                       *
1791
c     ****************************************************************
1792
 
1793
      implicit none
1794
 
1795
      integer mxnv
1796
 
1797
c     ----------------------------------------------------------------
1798
c     The statement below must match that in subroutine DefSPolyBndry
1799
c     ----------------------------------------------------------------
1800
 
1801
      parameter (mxnv=500)
1802
 
1803
      real tlonv(mxnv),vlat_c(mxnv),vlon_c(mxnv),xlat_c,xlon_c
1804
      real plat,plon,vAlat,vAlon,vBlat,vBlon,tlonA,tlonB,tlonP
1805
      real tlon_X,tlon_P,tlon_B,dellon
1806
      integer i,ibndry,nv_c,location,icross,ibrngAB,ibrngAP,ibrngPB
1807
      integer ibrng_BX,ibrng_BP,istrike
1808
 
1809
      common/spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry
1810
 
1811
      if (ibndry.eq.0) then     ! user has never defined the bndry
1812
         print*,'Subroutine LctPtRelBndry detects user error:'
1813
         print*,'Subroutine DefSPolyBndry must be called before'
1814
         print*,'subroutine LctPtRelBndry can be called'
1815
         stop
1816
      endif
1817
 
1818
      if (plat.eq.(-xlat_c)) then
1819
         dellon=plon-xlon_c
1820
         if (dellon.lt.(-180.)) dellon=dellon+360.
1821
         if (dellon.gt.+180.) dellon=dellon-360.
1822
         if ((dellon.eq.+180.0).or.(dellon.eq.-180.)) then
1823
            print*,'Warning: LctPtRelBndry detects case P antipodal
1824
     >           to X'
1825
            print*,'location of P relative to S is undetermined'
1826
            location=3
1827
            return
1828
         endif
1829
      endif 
1830
 
1831
      location=0                ! default ( P is outside S)
1832
      icross=0                  ! initialize counter
1833
 
1834
      if ((plat.eq.xlat_c).and.(plon.eq.xlon_c)) then
1835
         location=1
1836
         return
1837
      endif
1838
 
1839
 
1840
      call TrnsfmLon (xlat_c,xlon_c,plat,plon,tlonP)
1841
 
1842
      do i=1,nv_c              ! start of loop over sides of S 
1843
 
1844
         vAlat=vlat_c(i)
1845
         vAlon=vlon_c(i)
1846
         tlonA=tlonv(i)
1847
 
1848
         if (i.lt.nv_c) then
1849
            vBlat=vlat_c(i+1)
1850
            vBlon=vlon_c(i+1)
1851
            tlonB=tlonv(i+1)
1852
         else
1853
            vBlat=vlat_c(1)
1854
            vBlon=vlon_c(1)
1855
            tlonB=tlonv(1)
1856
         endif
1857
 
1858
         istrike=0
1859
 
1860
         if (tlonP.eq.tlonA) then
1861
            istrike=1
1862
         else
1863
            call EastOrWest(tlonA,tlonB,ibrngAB)
1864
            call EastOrWest(tlonA,tlonP,ibrngAP)
1865
            call EastOrWest(tlonP,tlonB,ibrngPB)
1866
 
1867
 
1868
            if((ibrngAP.eq.ibrngAB).and.(ibrngPB.eq.ibrngAB)) istrike=1
1869
         endif
1870
 
1871
 
1872
         if (istrike.eq.1) then
1873
 
1874
            if ((plat.eq.vAlat).and.(plon.eq.vAlon)) then
1875
               location=2       ! P lies on a vertex of S
1876
               return
1877
            endif
1878
            call TrnsfmLon(vAlat,vAlon,xlat_c,xlon_c,tlon_X)
1879
            call TrnsfmLon(vAlat,vAlon,vBlat,vBlon,tlon_B)
1880
            call TrnsfmLon(vAlat,vAlon,plat,plon,tlon_P)
1881
 
1882
            if (tlon_P.eq.tlon_B) then
1883
               location=2       ! P lies on side of S
1884
               return 
1885
            else
1886
               call EastOrWest(tlon_B,tlon_X,ibrng_BX)
1887
               call EastOrWest(tlon_B,tlon_P,ibrng_BP)
1888
               if(ibrng_BX.eq.(-ibrng_BP)) icross=icross+1
1889
            endif
1890
 
1891
         endif
1892
      enddo                     ! end of loop over the sides of S
1893
 
1894
 
1895
c     if the arc XP crosses the boundary S an even number of times then P
1896
c     is in S
1897
 
1898
      if (mod(icross,2).eq.0) location=1
1899
 
1900
      return
1901
 
1902
      end
1903
 
1904
 
1905
c     ****************************************************************
1906
 
1907
      subroutine TrnsfmLon(plat,plon,qlat,qlon,tranlon)
1908
 
1909
c     ****************************************************************
1910
c     * This subroutine is required by subroutines DefSPolyBndry &   *
1911
c     * LctPtRelBndry. It finds the 'longitude' of point Q in a      *
1912
c     * geographic coordinate system for which point P acts as a     *
1913
c     * 'north pole'. SENT: plat,plon,qlat,qlon, in degrees.         *
1914
c     * RETURNED: tranlon, in degrees.                               *
1915
c     ****************************************************************
1916
 
1917
      implicit none
1918
 
1919
      real pi,dtr,plat,plon,qlat,qlon,tranlon,t,b
1920
      parameter (pi=3.141592654,dtr=pi/180.0)
1921
 
1922
      if (plat.eq.90.) then
1923
         tranlon=qlon
1924
      else
1925
         t=sin((qlon-plon)*dtr)*cos(qlat*dtr)
1926
         b=sin(dtr*qlat)*cos(plat*dtr)-cos(qlat*dtr)*sin(plat*dtr)
1927
     >    *cos((qlon-plon)*dtr)
1928
         tranlon=atan2(t,b)/dtr
1929
      endif
1930
 
1931
      return
1932
      end
1933
 
1934
c     ****************************************************************
1935
 
1936
      subroutine EastOrWest(clon,dlon,ibrng)
1937
 
1938
c     ****************************************************************
1939
c     * This subroutine is required by subroutine LctPtRelBndry.     *
1940
c     * This routine determines if in travelling the shortest path   *
1941
c     * from point C (at longitude clon) to point D (at longitude    *
1942
c     * dlon) one is heading east, west or neither.                  *
1943
c     * SENT: clon,dlon; in degrees. RETURNED: ibrng                 *
1944
c     * (1=east,-1=west, 0=neither).                                 *
1945
c     ****************************************************************
1946
 
1947
      implicit none
1948
      real clon,dlon,del
1949
      integer ibrng
1950
      del=dlon-clon
1951
      if (del.gt.180.) del=del-360.
1952
      if (del.lt.-180.) del=del+360.
1953
      if ((del.gt.0.0).and.(del.ne.180.)) then
1954
         ibrng=-1               ! (D is west of C)
1955
      elseif ((del.lt.0.0).and.(del.ne.-180.)) then
1956
         ibrng=+1               ! (D is east of C)
1957
      else
1958
         ibrng=0                ! (D north or south of C)
1959
      endif
1960
      return
1961
      end