Subversion Repositories lagranto.um

Rev

Go to most recent revision | Details | Last modification | View Log | RSS feed

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