| 3 |
michaesp |
1 |
|
|
|
2 |
c -----------------------------------------------------------------------
|
|
|
3 |
c Nearest distance to a lat/lon point
|
|
|
4 |
c -----------------------------------------------------------------------
|
|
|
5 |
|
|
|
6 |
subroutine get_proxy (proxy,clon,clat,tra,ntra,ntim,ncol,radius)
|
|
|
7 |
|
|
|
8 |
c Given a set of trajectories tra(ntra,ntim,ncol), find the closest point
|
|
|
9 |
c to clon/clat and return the footprint of this trajectory at this
|
|
|
10 |
c point: proxy(ntra,ncol+1); the parameter <ncol+1> of the trajectory will
|
|
|
11 |
c become the minimum distance.
|
|
|
12 |
|
|
|
13 |
implicit none
|
|
|
14 |
|
|
|
15 |
c Declaration of parameters
|
|
|
16 |
integer ntra,ntim,ncol
|
|
|
17 |
real tra(ntra,ntim,ncol)
|
|
|
18 |
real proxy(ntra,ncol+1)
|
|
|
19 |
real radius
|
|
|
20 |
real clon,clat
|
|
|
21 |
|
|
|
22 |
c Flag for tests
|
|
|
23 |
integer test
|
|
|
24 |
parameter (test=0)
|
|
|
25 |
character*80 testfile
|
|
|
26 |
parameter (testfile='TEST.nc')
|
|
|
27 |
|
|
|
28 |
c Number of grid points for the radius mask
|
|
|
29 |
integer nmax
|
|
|
30 |
parameter (nmax = 400)
|
|
|
31 |
real degkm
|
|
|
32 |
parameter (degkm = 111.)
|
|
|
33 |
|
|
|
34 |
c Number of binary search steps
|
|
|
35 |
integer nbin
|
|
|
36 |
parameter (nbin=10)
|
|
|
37 |
|
|
|
38 |
c Auxiliry variables
|
|
|
39 |
integer i,j,k
|
|
|
40 |
real lon,lat,rlon,rlat
|
|
|
41 |
real dist
|
|
|
42 |
character*80 varname
|
|
|
43 |
integer rnx,rny
|
|
|
44 |
real rmask (nmax,nmax)
|
|
|
45 |
real rcount(nmax,nmax)
|
|
|
46 |
real rxmin,rymin,rdx,rdy
|
|
|
47 |
integer indx,indy
|
|
|
48 |
integer isnear
|
|
|
49 |
real near,near0,near1,nearm
|
|
|
50 |
integer j0,j1
|
|
|
51 |
real r0,r1,rm,frac
|
|
|
52 |
real rlon0,rlon1,rlonm,rlat0,rlat1,rlatm
|
|
|
53 |
|
|
|
54 |
c Externals
|
|
|
55 |
real sdis
|
|
|
56 |
external sdis
|
|
|
57 |
|
|
|
58 |
c Transform into clon/clat centered system
|
|
|
59 |
do i=1,ntra
|
|
|
60 |
do j=1,ntim
|
|
|
61 |
lon = tra(i,j,2)
|
|
|
62 |
lat = tra(i,j,3)
|
|
|
63 |
call getenvir_f (clon,clat,0.,rlon,rlat,lon,lat,1)
|
|
|
64 |
tra(i,j,2) = rlon
|
|
|
65 |
tra(i,j,3) = rlat
|
|
|
66 |
enddo
|
|
|
67 |
enddo
|
|
|
68 |
|
|
|
69 |
c Init the radius mask - define the grid resolution
|
|
|
70 |
rdx = 2.5 * radius / degkm / real(nmax)
|
|
|
71 |
rdy = 2.5 * radius / degkm / real(nmax)
|
|
|
72 |
rnx = nmax
|
|
|
73 |
rny = nmax
|
|
|
74 |
rxmin = -real(rnx/2)*rdx
|
|
|
75 |
rymin = -real(rny/2)*rdy
|
|
|
76 |
do i=1,rnx
|
|
|
77 |
do j=1,rny
|
|
|
78 |
rlon = rxmin + real(i-1) * rdx
|
|
|
79 |
rlat = rymin + real(j-1) * rdy
|
|
|
80 |
dist = sdis(rlon,rlat,0.,0.)
|
|
|
81 |
if ( dist.lt.radius ) then
|
|
|
82 |
rmask(i,j) = dist
|
|
|
83 |
else
|
|
|
84 |
rmask(i,j) = radius
|
|
|
85 |
endif
|
|
|
86 |
enddo
|
|
|
87 |
enddo
|
|
|
88 |
|
|
|
89 |
c Init counter for test
|
|
|
90 |
if ( test.eq.1 ) then
|
|
|
91 |
do i=1,rnx
|
|
|
92 |
do j=1,rny
|
|
|
93 |
rcount(i,j) = 0.
|
|
|
94 |
enddo
|
|
|
95 |
enddo
|
|
|
96 |
endif
|
|
|
97 |
|
|
|
98 |
|
|
|
99 |
c Loop over all trajectories
|
|
|
100 |
do i=1,ntra
|
|
|
101 |
|
|
|
102 |
c Decide whether trajectory comes close to point
|
|
|
103 |
isnear = 0
|
|
|
104 |
near = radius
|
|
|
105 |
do j=1,ntim
|
|
|
106 |
indx = nint( ( tra(i,j,2) - rxmin ) / rdx + 1. )
|
|
|
107 |
indy = nint( ( tra(i,j,3) - rymin ) / rdy + 1. )
|
|
|
108 |
if ( (indx.ge.1).and.(indx.le.rnx).and.
|
|
|
109 |
> (indy.ge.1).and.(indy.le.rny) )
|
|
|
110 |
> then
|
|
|
111 |
if (rmask(indx,indy).lt.near) then
|
|
|
112 |
near = rmask(indx,indy)
|
|
|
113 |
isnear = j
|
|
|
114 |
endif
|
|
|
115 |
endif
|
|
|
116 |
enddo
|
|
|
117 |
|
|
|
118 |
c No close point was found - go to next trajectory
|
|
|
119 |
do k=1,ncol
|
|
|
120 |
proxy(i,k) = tra(i,1,k)
|
|
|
121 |
enddo
|
|
|
122 |
proxy(i,ncol+1) = radius
|
|
|
123 |
if ( isnear.eq.0 ) goto 310
|
|
|
124 |
|
|
|
125 |
c Get the exact position and time with binary search
|
|
|
126 |
j0 = isnear - 1
|
|
|
127 |
if ( j0.eq.0 ) j0 = 1
|
|
|
128 |
rlon0 = tra(i,j0,2)
|
|
|
129 |
rlat0 = tra(i,j0,3)
|
|
|
130 |
near0 = sdis(rlon0,rlat0,0.,0.)
|
|
|
131 |
r0 = real(j0)
|
|
|
132 |
|
|
|
133 |
j1 = isnear + 1
|
|
|
134 |
if ( j1.gt.ntim ) j1 = ntim
|
|
|
135 |
rlon1 = tra(i,j1,2)
|
|
|
136 |
rlat1 = tra(i,j1,3)
|
|
|
137 |
near1 = sdis(tra(i,j1,2),tra(i,j1,3),0.,0.)
|
|
|
138 |
r1 = real(j1)
|
|
|
139 |
|
|
|
140 |
do k=1,nbin
|
|
|
141 |
rm = 0.5 * ( r0 + r1 )
|
|
|
142 |
rlatm = 0.5 * ( rlat0 + rlat1 )
|
|
|
143 |
rlonm = 0.5 * ( rlon0 + rlon1 )
|
|
|
144 |
nearm = sdis(rlonm,rlatm,0.,0.)
|
|
|
145 |
if ( near0.lt.near1 ) then
|
|
|
146 |
r1 = rm
|
|
|
147 |
rlon1 = rlonm
|
|
|
148 |
rlat1 = rlatm
|
|
|
149 |
near1 = nearm
|
|
|
150 |
else
|
|
|
151 |
r0 = rm
|
|
|
152 |
rlon0 = rlonm
|
|
|
153 |
rlat0 = rlatm
|
|
|
154 |
near0 = nearm
|
|
|
155 |
endif
|
|
|
156 |
enddo
|
|
|
157 |
|
|
|
158 |
c Now get the final distance and position
|
|
|
159 |
proxy(i,ncol+1) = nearm
|
|
|
160 |
if ( test.eq.1 ) then
|
|
|
161 |
indx = nint( ( rlonm - rxmin ) / rdx + 1. )
|
|
|
162 |
indy = nint( ( rlatm - rymin ) / rdy + 1. )
|
|
|
163 |
if ( (indx.ge.1).and.(indx.le.rnx).and.
|
|
|
164 |
> (indy.ge.1).and.(indy.le.rny) )
|
|
|
165 |
> then
|
|
|
166 |
rcount(indx,indy) = rcount(indx,indy) + 1.
|
|
|
167 |
endif
|
|
|
168 |
endif
|
|
|
169 |
|
|
|
170 |
c Get all values at exactly this position
|
|
|
171 |
j0 = int(rm)
|
|
|
172 |
frac = rm - real(j0)
|
|
|
173 |
j1 = j0 + 1
|
|
|
174 |
if (j1.gt.ntim ) j1 = ntim
|
|
|
175 |
do k=1,ncol
|
|
|
176 |
proxy(i,k) = (1.-frac)*tra(i,j0,k)+frac*tra(i,j1,k)
|
|
|
177 |
enddo
|
|
|
178 |
|
|
|
179 |
c Next trajectory
|
|
|
180 |
310 continue
|
|
|
181 |
|
|
|
182 |
enddo
|
|
|
183 |
|
|
|
184 |
c Transform into geo system
|
|
|
185 |
do i=1,ntra
|
|
|
186 |
rlon = proxy(i,2)
|
|
|
187 |
rlat = proxy(i,3)
|
|
|
188 |
call getenvir_f (clon,clat,0.,lon,lat,rlon,rlat,1)
|
|
|
189 |
proxy(i,2) = lon
|
|
|
190 |
proxy(i,3) = lat
|
|
|
191 |
enddo
|
|
|
192 |
|
|
|
193 |
c Write radius mask to netCDF file for tests
|
|
|
194 |
if ( test.eq.1 ) then
|
|
|
195 |
varname = 'MASK'
|
|
|
196 |
call writecdf2D(testfile,varname,rmask,0.,
|
|
|
197 |
> rdx,rdy,rxmin,rymin,rnx,rny,1,1)
|
|
|
198 |
varname = 'COUNT'
|
|
|
199 |
call writecdf2D(testfile,varname,rcount,0.,
|
|
|
200 |
> rdx,rdy,rxmin,rymin,rnx,rny,0,1)
|
|
|
201 |
endif
|
|
|
202 |
|
|
|
203 |
end
|
|
|
204 |
|
|
|
205 |
|
|
|
206 |
|
|
|
207 |
c --------------------------------------------------------------------------------
|
|
|
208 |
c Forward coordinate transformation (True lon/lat -> Rotated lon/lat)
|
|
|
209 |
c --------------------------------------------------------------------------------
|
|
|
210 |
|
|
|
211 |
SUBROUTINE getenvir_f (clon,clat,rotation,
|
|
|
212 |
> rlon,rlat,lon,lat,n)
|
|
|
213 |
|
|
|
214 |
implicit none
|
|
|
215 |
|
|
|
216 |
c Declaration of input and output parameters
|
|
|
217 |
integer n
|
|
|
218 |
real clon,clat,rotation
|
|
|
219 |
real lon(n), lat(n)
|
|
|
220 |
real rlon(n),rlat(n)
|
|
|
221 |
|
|
|
222 |
c Auxiliary variables
|
|
|
223 |
real pollon,pollat
|
|
|
224 |
integer i
|
|
|
225 |
real rlon1(n),rlat1(n)
|
|
|
226 |
|
|
|
227 |
c Externals
|
|
|
228 |
real lmtolms,phtophs
|
|
|
229 |
external lmtolms,phtophs
|
|
|
230 |
|
|
|
231 |
c First rotation
|
|
|
232 |
pollon=clon-180.
|
|
|
233 |
if (pollon.lt.-180.) pollon=pollon+360.
|
|
|
234 |
pollat=90.-clat
|
|
|
235 |
do i=1,n
|
|
|
236 |
|
|
|
237 |
c First rotation
|
|
|
238 |
pollon=clon-180.
|
|
|
239 |
if (pollon.lt.-180.) pollon=pollon+360.
|
|
|
240 |
pollat=90.-clat
|
|
|
241 |
rlon1(i)=lmtolms(lat(i),lon(i),pollat,pollon)
|
|
|
242 |
rlat1(i)=phtophs(lat(i),lon(i),pollat,pollon)
|
|
|
243 |
|
|
|
244 |
c Second coordinate transformation
|
|
|
245 |
pollon=-180.
|
|
|
246 |
pollat=90.+rotation
|
|
|
247 |
rlon(i)=90.+lmtolms(rlat1(i),rlon1(i)-90.,pollat,pollon)
|
|
|
248 |
rlat(i)=phtophs(rlat1(i),rlon1(i)-90.,pollat,pollon)
|
|
|
249 |
|
|
|
250 |
enddo
|
|
|
251 |
|
|
|
252 |
END
|
|
|
253 |
|
|
|
254 |
c --------------------------------------------------------------------------------
|
|
|
255 |
c Backward coordinate transformation (Rotated lon/lat -> True lon/lat)
|
|
|
256 |
c --------------------------------------------------------------------------------
|
|
|
257 |
|
|
|
258 |
SUBROUTINE getenvir_b (clon,clat,rotation,
|
|
|
259 |
> lon,lat,rlon,rlat,n)
|
|
|
260 |
|
|
|
261 |
implicit none
|
|
|
262 |
|
|
|
263 |
c Declaration of input and output parameters
|
|
|
264 |
integer n
|
|
|
265 |
real clon,clat,rotation
|
|
|
266 |
real lon(n), lat(n)
|
|
|
267 |
real rlon(n),rlat(n)
|
|
|
268 |
|
|
|
269 |
c Auxiliary variables
|
|
|
270 |
real pollon,pollat
|
|
|
271 |
integer i
|
|
|
272 |
real rlon1(n),rlat1(n)
|
|
|
273 |
|
|
|
274 |
c Externals
|
|
|
275 |
real lmstolm,phstoph
|
|
|
276 |
external lmstolm,phstoph
|
|
|
277 |
|
|
|
278 |
c First coordinate transformation (make the local coordinate system parallel to equator)
|
|
|
279 |
pollon=-180.
|
|
|
280 |
pollat=90.+rotation
|
|
|
281 |
do i=1,n
|
|
|
282 |
rlon1(i)=90.+lmstolm(rlat(i),rlon(i)-90.,pollat,pollon)
|
|
|
283 |
rlat1(i)=phstoph(rlat(i),rlon(i)-90.,pollat,pollon)
|
|
|
284 |
enddo
|
|
|
285 |
|
|
|
286 |
c Second coordinate transformation (make the local coordinate system parallel to equator)
|
|
|
287 |
pollon=clon-180.
|
|
|
288 |
if (pollon.lt.-180.) pollon=pollon+360.
|
|
|
289 |
pollat=90.-clat
|
|
|
290 |
do i=1,n
|
|
|
291 |
lon(i)=lmstolm(rlat1(i),rlon1(i),pollat,pollon)
|
|
|
292 |
lat(i)=phstoph(rlat1(i),rlon1(i),pollat,pollon)
|
|
|
293 |
enddo
|
|
|
294 |
|
|
|
295 |
END
|
|
|
296 |
|
|
|
297 |
c -----------------------------------------------------------------------
|
|
|
298 |
c Spherical distance between lat/lon points
|
|
|
299 |
c -----------------------------------------------------------------------
|
|
|
300 |
|
|
|
301 |
real function sdis(xp,yp,xq,yq)
|
|
|
302 |
c
|
|
|
303 |
c calculates spherical distance (in km) between two points given
|
|
|
304 |
c by their spherical coordinates (xp,yp) and (xq,yq), respectively.
|
|
|
305 |
c
|
|
|
306 |
real re
|
|
|
307 |
parameter (re=6370.)
|
|
|
308 |
real pi180
|
|
|
309 |
parameter (pi180=3.14159/180.)
|
|
|
310 |
real xp,yp,xq,yq,arg
|
|
|
311 |
|
|
|
312 |
arg=sin(pi180*yp)*sin(pi180*yq)+
|
|
|
313 |
> cos(pi180*yp)*cos(pi180*yq)*cos(pi180*(xp-xq))
|
|
|
314 |
if (arg.lt.-1.) arg=-1.
|
|
|
315 |
if (arg.gt.1.) arg=1.
|
|
|
316 |
|
|
|
317 |
sdis=re*acos(arg)
|
|
|
318 |
|
|
|
319 |
end
|
|
|
320 |
|