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