1 | !*********************************************************************** |
---|
2 | !* Copyright 2012,2013 * |
---|
3 | !* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine, * |
---|
4 | !* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,* |
---|
5 | !* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso, * |
---|
6 | !* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann, * |
---|
7 | !* * |
---|
8 | !* This file is part of FLEXPART WRF * |
---|
9 | !* * |
---|
10 | !* FLEXPART is free software: you can redistribute it and/or modify * |
---|
11 | !* it under the terms of the GNU General Public License as published by* |
---|
12 | !* the Free Software Foundation, either version 3 of the License, or * |
---|
13 | !* (at your option) any later version. * |
---|
14 | !* * |
---|
15 | !* FLEXPART is distributed in the hope that it will be useful, * |
---|
16 | !* but WITHOUT ANY WARRANTY; without even the implied warranty of * |
---|
17 | !* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * |
---|
18 | !* GNU General Public License for more details. * |
---|
19 | !* * |
---|
20 | !* You should have received a copy of the GNU General Public License * |
---|
21 | !* along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * |
---|
22 | !*********************************************************************** |
---|
23 | !********************************************************************** |
---|
24 | ! FLEXPART_WRF SOURCE FILE MAP_PROJ_WRF - CONTAINS * |
---|
25 | ! subroutine xyindex_to_ll_wrf * |
---|
26 | ! subroutine xymeter_to_ll_wrf * |
---|
27 | ! subroutine ll_to_xyindex_wrf * |
---|
28 | ! subroutine ll_to_xymeter_wrf * |
---|
29 | ! subroutine test_xyindex_to_ll_wrf * |
---|
30 | ! * |
---|
31 | !********************************************************************** |
---|
32 | ! * |
---|
33 | ! AUTHOR: R. Easter & J. Fast, PNNL * |
---|
34 | ! DATE: Dec 2005 * |
---|
35 | ! LAST UPDATE: same * |
---|
36 | ! * |
---|
37 | !********************************************************************** |
---|
38 | ! * |
---|
39 | ! DESCRIPTION: * |
---|
40 | ! * |
---|
41 | ! Converts between "grid index" or "grid meter" coordinates * |
---|
42 | ! and latitude-longitude (degrees) * |
---|
43 | ! * |
---|
44 | !********************************************************************** |
---|
45 | |
---|
46 | |
---|
47 | !----------------------------------------------------------------------- |
---|
48 | subroutine xyindex_to_ll_wrf( lgrid, x_i, y_j, x_lon, y_lat ) |
---|
49 | ! |
---|
50 | ! calculates longitude/latitude from xy "grid index" coordinates |
---|
51 | ! |
---|
52 | ! arguments |
---|
53 | ! lgrid - input - grid identifier |
---|
54 | ! x_i, y_j - input = "grid index" coordinates on grid "lgrid". |
---|
55 | ! x_i ranges from 0 to nx-1. y_j ranges from 0 to nj-1. |
---|
56 | ! x_lon, y_lat - output = longitude and latitude in degrees |
---|
57 | ! |
---|
58 | ! *** note *** |
---|
59 | ! if x_i is outside [-grace, +grace+nx] |
---|
60 | ! or y_j is outside [-grace, +grace+ny] |
---|
61 | ! the routine writes an error message and halts |
---|
62 | ! |
---|
63 | ! include 'includepar' |
---|
64 | ! include 'includecom' |
---|
65 | use par_mod |
---|
66 | use com_mod |
---|
67 | |
---|
68 | ! arguments |
---|
69 | implicit none |
---|
70 | integer :: lgrid |
---|
71 | real :: x_i, y_j, x_lon, y_lat |
---|
72 | |
---|
73 | ! local variables |
---|
74 | integer :: ia, ib, ja, jb |
---|
75 | real :: dumxi, dumyj |
---|
76 | real :: fx, fy |
---|
77 | real,parameter :: grace=2.01 |
---|
78 | |
---|
79 | |
---|
80 | ! first check that x_i & y_j are "within bounds" |
---|
81 | if (lgrid .le. 0) then |
---|
82 | if ((x_i .lt. -grace) .or. (x_i .gt. grace+real(nx-1)) .or. & |
---|
83 | (y_j .lt. -grace) .or. (y_j .gt. grace+real(ny-1))) then |
---|
84 | write(*,'(/a/a,i4,1p,2e12.3)') & |
---|
85 | '*** xyindex_to_ll_wrf error -- bad inputs', & |
---|
86 | ' lgrid, x_i, y_j =', lgrid, x_i, y_j |
---|
87 | ! stop |
---|
88 | end if |
---|
89 | else |
---|
90 | if ((x_i .lt. -grace ) .or. & |
---|
91 | (x_i .gt. grace+real(nxn(lgrid)-1)) .or. & |
---|
92 | (y_j .lt. -grace ) .or. & |
---|
93 | (y_j .gt. grace+real(nyn(lgrid)-1))) then |
---|
94 | write(*,'(/a/a,i4,1p,2e12.3)') & |
---|
95 | '*** xyindex_to_ll_wrf error -- bad inputs', & |
---|
96 | ' lgrid, x_i, y_j =', lgrid, x_i, y_j |
---|
97 | ! stop |
---|
98 | end if |
---|
99 | end if |
---|
100 | |
---|
101 | if (map_proj_method .gt. 0) then |
---|
102 | if (map_proj_id .eq. 1) goto 2000 |
---|
103 | if (map_proj_id .eq. 2) goto 2000 |
---|
104 | if (map_proj_id .eq. 3) goto 2000 |
---|
105 | if (map_proj_id .eq. 4) goto 2000 |
---|
106 | end if |
---|
107 | |
---|
108 | if (lgrid .le. 0) then |
---|
109 | goto 5000 |
---|
110 | else |
---|
111 | goto 6000 |
---|
112 | end if |
---|
113 | |
---|
114 | ! |
---|
115 | ! use map projection routines in map_proj_wrf_subaa.f |
---|
116 | ! |
---|
117 | 2000 continue |
---|
118 | if (lgrid .le. 0) then |
---|
119 | dumxi = 1.0 + x_i |
---|
120 | dumyj = 1.0 + y_j |
---|
121 | else |
---|
122 | dumxi = 1.0 + (x_i/xresoln(lgrid)) + xln(lgrid) |
---|
123 | dumyj = 1.0 + (y_j/yresoln(lgrid)) + yln(lgrid) |
---|
124 | end if |
---|
125 | call ij_to_latlon( dumxi, dumyj, y_lat, x_lon ) |
---|
126 | return |
---|
127 | |
---|
128 | |
---|
129 | ! |
---|
130 | ! do interpolation using the outer grid xlon2d,ylat2d |
---|
131 | ! |
---|
132 | 5000 continue |
---|
133 | if (x_i .le. 0.0) then |
---|
134 | ia = 0 |
---|
135 | else if (x_i .ge. real(nxmin1)) then |
---|
136 | ia = nxmin1 |
---|
137 | else |
---|
138 | ia = ifix( x_i ) |
---|
139 | end if |
---|
140 | fx = x_i - ia |
---|
141 | fx = max( -2.0, min( fx, 3.0 ) ) |
---|
142 | ib = ia + 1 |
---|
143 | |
---|
144 | if (y_j .le. 0.0) then |
---|
145 | ja = 0 |
---|
146 | else if (y_j .ge. real(nymin1)) then |
---|
147 | ja = nymin1 |
---|
148 | else |
---|
149 | ja = ifix( y_j ) |
---|
150 | end if |
---|
151 | fy = y_j - ja |
---|
152 | fy = max( -2.0, min( fy, 3.0 ) ) |
---|
153 | jb = ja + 1 |
---|
154 | |
---|
155 | x_lon = xlon2d(ia,ja)*(1.0-fx)*(1.0-fy) + & |
---|
156 | xlon2d(ia,jb)*(1.0-fx)*fy + & |
---|
157 | xlon2d(ib,ja)*fx*(1.0-fy) + & |
---|
158 | xlon2d(ib,jb)*fx*fy |
---|
159 | |
---|
160 | |
---|
161 | y_lat = ylat2d(ia,ja)*(1.0-fx)*(1.0-fy) + & |
---|
162 | ylat2d(ia,jb)*(1.0-fx)*fy + & |
---|
163 | ylat2d(ib,ja)*fx*(1.0-fy) + & |
---|
164 | ylat2d(ib,jb)*fx*fy |
---|
165 | |
---|
166 | return |
---|
167 | |
---|
168 | |
---|
169 | ! |
---|
170 | ! do interpolation using the nested grid xlon2dn,ylat2dn |
---|
171 | ! |
---|
172 | 6000 continue |
---|
173 | if (x_i .le. 0.0) then |
---|
174 | ia = 0 |
---|
175 | else if (x_i .ge. real(nxn(lgrid)-1)) then |
---|
176 | ia = nxn(lgrid)-1 |
---|
177 | else |
---|
178 | ia = ifix( x_i ) |
---|
179 | end if |
---|
180 | fx = x_i - ia |
---|
181 | fx = max( -2.0, min( fx, 3.0 ) ) |
---|
182 | ib = ia + 1 |
---|
183 | |
---|
184 | if (y_j .le. 0.0) then |
---|
185 | ja = 0 |
---|
186 | else if (y_j .ge. real(nyn(lgrid)-1)) then |
---|
187 | ja = nyn(lgrid)-1 |
---|
188 | else |
---|
189 | ja = ifix( y_j ) |
---|
190 | end if |
---|
191 | fy = y_j - ja |
---|
192 | fy = max( -2.0, min( fy, 3.0 ) ) |
---|
193 | jb = ja + 1 |
---|
194 | |
---|
195 | x_lon = xlon2dn(ia,ja,lgrid)*(1.0-fx)*(1.0-fy) + & |
---|
196 | xlon2dn(ia,jb,lgrid)*(1.0-fx)*fy + & |
---|
197 | xlon2dn(ib,ja,lgrid)*fx*(1.0-fy) + & |
---|
198 | xlon2dn(ib,jb,lgrid)*fx*fy |
---|
199 | |
---|
200 | |
---|
201 | y_lat = ylat2dn(ia,ja,lgrid)*(1.0-fx)*(1.0-fy) + & |
---|
202 | ylat2dn(ia,jb,lgrid)*(1.0-fx)*fy + & |
---|
203 | ylat2dn(ib,ja,lgrid)*fx*(1.0-fy) + & |
---|
204 | ylat2dn(ib,jb,lgrid)*fx*fy |
---|
205 | |
---|
206 | return |
---|
207 | |
---|
208 | end subroutine xyindex_to_ll_wrf |
---|
209 | |
---|
210 | |
---|
211 | !----------------------------------------------------------------------- |
---|
212 | subroutine xyindex_to_ll_wrf_out(lgrid,x_i,y_j,x_lon,y_lat) |
---|
213 | ! |
---|
214 | ! calculates longitude/latitude from xy "grid index" coordinates |
---|
215 | ! |
---|
216 | ! arguments |
---|
217 | ! lgrid - input - grid identifier |
---|
218 | ! x_i, y_j - input = "grid index" coordinates on grid "lgrid". |
---|
219 | ! x_i ranges from 0 to nx-1. y_j ranges from 0 to nj-1. |
---|
220 | ! x_lon, y_lat - output = longitude and latitude in degrees |
---|
221 | ! |
---|
222 | ! *** note *** |
---|
223 | ! if x_i is outside [-grace, +grace+nx] |
---|
224 | ! or y_j is outside [-grace, +grace+ny] |
---|
225 | ! the routine writes an error message and halts |
---|
226 | ! |
---|
227 | ! include 'includepar' |
---|
228 | ! include 'includecom' |
---|
229 | use par_mod |
---|
230 | use com_mod |
---|
231 | |
---|
232 | ! arguments |
---|
233 | implicit none |
---|
234 | integer :: lgrid |
---|
235 | real :: x_i, y_j, x_lon, y_lat |
---|
236 | |
---|
237 | ! local variables |
---|
238 | integer :: ia, ib, ja, jb |
---|
239 | real :: dumxi, dumyj |
---|
240 | real :: fx, fy |
---|
241 | real,parameter :: grace=2.01 |
---|
242 | |
---|
243 | |
---|
244 | ! first check that x_i & y_j are "within bounds" |
---|
245 | ! if (lgrid .le. 0) then |
---|
246 | ! if ((x_i .lt. -grace) .or. (x_i .gt. grace+real(nx-1)) .or. |
---|
247 | ! & (y_j .lt. -grace) .or. (y_j .gt. grace+real(ny-1))) then |
---|
248 | ! write(*,'(/a/a,i4,1p,2e12.3)') |
---|
249 | ! & '*** xyindex_to_ll_wrf error -- bad inputs', |
---|
250 | ! & ' lgrid, x_i, y_j =', lgrid, x_i, y_j |
---|
251 | ! stop |
---|
252 | ! end if |
---|
253 | ! else |
---|
254 | ! if ((x_i .lt. -grace ) .or. |
---|
255 | ! & (x_i .gt. grace+real(nxn(lgrid)-1)) .or. |
---|
256 | ! & (y_j .lt. -grace ) .or. |
---|
257 | ! & (y_j .gt. grace+real(nyn(lgrid)-1))) then |
---|
258 | ! write(*,'(/a/a,i4,1p,2e12.3)') |
---|
259 | ! & '*** xyindex_to_ll_wrf error -- bad inputs', |
---|
260 | ! & ' lgrid, x_i, y_j =', lgrid, x_i, y_j |
---|
261 | ! stop |
---|
262 | ! end if |
---|
263 | ! end if |
---|
264 | |
---|
265 | if (map_proj_method .gt. 0) then |
---|
266 | if (map_proj_id .eq. 1) goto 2000 |
---|
267 | if (map_proj_id .eq. 2) goto 2000 |
---|
268 | if (map_proj_id .eq. 3) goto 2000 |
---|
269 | if (map_proj_id .eq. 4) goto 2000 |
---|
270 | end if |
---|
271 | |
---|
272 | if (lgrid .le. 0) then |
---|
273 | goto 5000 |
---|
274 | else |
---|
275 | goto 6000 |
---|
276 | end if |
---|
277 | |
---|
278 | ! |
---|
279 | ! use map projection routines in map_proj_wrf_subaa.f |
---|
280 | ! |
---|
281 | 2000 continue |
---|
282 | if (lgrid .le. 0) then |
---|
283 | dumxi = 1.0 + x_i |
---|
284 | dumyj = 1.0 + y_j |
---|
285 | else |
---|
286 | dumxi = 1.0 + (x_i/xresoln(lgrid)) + xln(lgrid) |
---|
287 | dumyj = 1.0 + (y_j/yresoln(lgrid)) + yln(lgrid) |
---|
288 | end if |
---|
289 | call ij_to_latlon( dumxi, dumyj, y_lat, x_lon ) |
---|
290 | return |
---|
291 | |
---|
292 | |
---|
293 | ! |
---|
294 | ! do interpolation using the outer grid xlon2d,ylat2d |
---|
295 | ! |
---|
296 | 5000 continue |
---|
297 | if (x_i .le. 0.0) then |
---|
298 | ia = 0 |
---|
299 | else if (x_i .ge. real(nxmin1)) then |
---|
300 | ia = nxmin1 |
---|
301 | else |
---|
302 | ia = ifix( x_i ) |
---|
303 | end if |
---|
304 | fx = x_i - ia |
---|
305 | fx = max( -2.0, min( fx, 3.0 ) ) |
---|
306 | ib = ia + 1 |
---|
307 | |
---|
308 | if (y_j .le. 0.0) then |
---|
309 | ja = 0 |
---|
310 | else if (y_j .ge. real(nymin1)) then |
---|
311 | ja = nymin1 |
---|
312 | else |
---|
313 | ja = ifix( y_j ) |
---|
314 | end if |
---|
315 | fy = y_j - ja |
---|
316 | fy = max( -2.0, min( fy, 3.0 ) ) |
---|
317 | jb = ja + 1 |
---|
318 | |
---|
319 | x_lon = xlon2d(ia,ja)*(1.0-fx)*(1.0-fy) + & |
---|
320 | xlon2d(ia,jb)*(1.0-fx)*fy + & |
---|
321 | xlon2d(ib,ja)*fx*(1.0-fy) + & |
---|
322 | xlon2d(ib,jb)*fx*fy |
---|
323 | |
---|
324 | |
---|
325 | y_lat = ylat2d(ia,ja)*(1.0-fx)*(1.0-fy) + & |
---|
326 | ylat2d(ia,jb)*(1.0-fx)*fy + & |
---|
327 | ylat2d(ib,ja)*fx*(1.0-fy) + & |
---|
328 | ylat2d(ib,jb)*fx*fy |
---|
329 | |
---|
330 | return |
---|
331 | |
---|
332 | |
---|
333 | ! |
---|
334 | ! do interpolation using the nested grid xlon2dn,ylat2dn |
---|
335 | ! |
---|
336 | 6000 continue |
---|
337 | if (x_i .le. 0.0) then |
---|
338 | ia = 0 |
---|
339 | else if (x_i .ge. real(nxn(lgrid)-1)) then |
---|
340 | ia = nxn(lgrid)-1 |
---|
341 | else |
---|
342 | ia = ifix( x_i ) |
---|
343 | end if |
---|
344 | fx = x_i - ia |
---|
345 | fx = max( -2.0, min( fx, 3.0 ) ) |
---|
346 | ib = ia + 1 |
---|
347 | |
---|
348 | if (y_j .le. 0.0) then |
---|
349 | ja = 0 |
---|
350 | else if (y_j .ge. real(nyn(lgrid)-1)) then |
---|
351 | ja = nyn(lgrid)-1 |
---|
352 | else |
---|
353 | ja = ifix( y_j ) |
---|
354 | end if |
---|
355 | fy = y_j - ja |
---|
356 | fy = max( -2.0, min( fy, 3.0 ) ) |
---|
357 | jb = ja + 1 |
---|
358 | |
---|
359 | x_lon = xlon2dn(ia,ja,lgrid)*(1.0-fx)*(1.0-fy) + & |
---|
360 | xlon2dn(ia,jb,lgrid)*(1.0-fx)*fy + & |
---|
361 | xlon2dn(ib,ja,lgrid)*fx*(1.0-fy) + & |
---|
362 | xlon2dn(ib,jb,lgrid)*fx*fy |
---|
363 | |
---|
364 | |
---|
365 | y_lat = ylat2dn(ia,ja,lgrid)*(1.0-fx)*(1.0-fy) + & |
---|
366 | ylat2dn(ia,jb,lgrid)*(1.0-fx)*fy + & |
---|
367 | ylat2dn(ib,ja,lgrid)*fx*(1.0-fy) + & |
---|
368 | ylat2dn(ib,jb,lgrid)*fx*fy |
---|
369 | |
---|
370 | return |
---|
371 | |
---|
372 | end subroutine xyindex_to_ll_wrf_out |
---|
373 | |
---|
374 | |
---|
375 | !----------------------------------------------------------------------- |
---|
376 | subroutine xymeter_to_ll_wrf( xmeter, ymeter, x_lon, y_lat ) |
---|
377 | ! |
---|
378 | ! calculates longitude/latitude from xy "grid meter" coordinates |
---|
379 | ! |
---|
380 | ! arguments |
---|
381 | ! xmeter, ymeter - input = "grid meter" coordinates on the mother grid. |
---|
382 | ! x_lon, y_lat - output = longitude and latitude in degrees |
---|
383 | ! |
---|
384 | ! include 'includepar' |
---|
385 | ! include 'includecom' |
---|
386 | use par_mod |
---|
387 | use com_mod |
---|
388 | |
---|
389 | ! arguments |
---|
390 | implicit none |
---|
391 | real :: xmeter, ymeter, x_lon, y_lat |
---|
392 | |
---|
393 | ! local variables |
---|
394 | real :: x_i, y_j |
---|
395 | |
---|
396 | x_i = (xmeter - xmet0)/dx |
---|
397 | y_j = (ymeter - ymet0)/dy |
---|
398 | call xyindex_to_ll_wrf( 0, x_i, y_j, x_lon, y_lat ) |
---|
399 | |
---|
400 | return |
---|
401 | end subroutine xymeter_to_ll_wrf |
---|
402 | |
---|
403 | subroutine xymeter_to_ll_wrf_out(xmeter,ymeter,x_lon,y_lat) |
---|
404 | ! |
---|
405 | ! calculates longitude/latitude from xy "grid meter" coordinates |
---|
406 | ! |
---|
407 | ! arguments |
---|
408 | ! xmeter, ymeter - input = "grid meter" coordinates on the mother grid. |
---|
409 | ! x_lon, y_lat - output = longitude and latitude in degrees |
---|
410 | ! |
---|
411 | ! include 'includepar' |
---|
412 | ! include 'includecom' |
---|
413 | use par_mod |
---|
414 | use com_mod |
---|
415 | |
---|
416 | ! arguments |
---|
417 | implicit none |
---|
418 | real :: xmeter, ymeter, x_lon, y_lat |
---|
419 | |
---|
420 | ! local variables |
---|
421 | real :: x_i, y_j |
---|
422 | |
---|
423 | x_i = (xmeter - xmet0)/dx |
---|
424 | y_j = (ymeter - ymet0)/dy |
---|
425 | call xyindex_to_ll_wrf_out(0,x_i, y_j,x_lon,y_lat ) |
---|
426 | |
---|
427 | return |
---|
428 | end subroutine xymeter_to_ll_wrf_out |
---|
429 | |
---|
430 | !----------------------------------------------------------------------- |
---|
431 | subroutine ll_to_xyindex_wrf( x_lon, y_lat, x_i, y_j ) |
---|
432 | ! |
---|
433 | ! calculates xy "grid index" coordinates from longitude/latitude |
---|
434 | ! |
---|
435 | ! arguments |
---|
436 | ! x_lon, y_lat - input - longitude and latitude in degrees |
---|
437 | ! x_i, y_j - output = "grid index" coordinates on the mother grid. |
---|
438 | ! x_i ranges from 0 to nx-1. y_j ranges from 0 to nj-1. |
---|
439 | ! |
---|
440 | ! include 'includepar' |
---|
441 | ! include 'includecom' |
---|
442 | use par_mod |
---|
443 | use com_mod |
---|
444 | implicit none |
---|
445 | |
---|
446 | ! arguments |
---|
447 | real :: x_i, y_j, x_lon, y_lat |
---|
448 | |
---|
449 | ! local variables |
---|
450 | integer :: i, ii, ia, ib, ip, ipass, iijjlohi |
---|
451 | integer :: j, jj, ja, jb, jp |
---|
452 | real :: dumcos, dumlat, dumlon |
---|
453 | real :: dumr2, dumr2min |
---|
454 | real :: dumxi, dumyj |
---|
455 | real :: flo, fhi |
---|
456 | real :: glo, ghi |
---|
457 | real,parameter :: grace=2.01 |
---|
458 | real :: xxcen, xxcenb, yycen, yycenb |
---|
459 | real :: xxyydel, xxyydelmin |
---|
460 | real :: x_lon_sv, y_lat_sv |
---|
461 | |
---|
462 | |
---|
463 | |
---|
464 | x_lon_sv = x_lon |
---|
465 | y_lat_sv = y_lat |
---|
466 | if (map_proj_method .gt. 0) then |
---|
467 | if (map_proj_id .eq. 1) goto 2000 |
---|
468 | if (map_proj_id .eq. 2) goto 2000 |
---|
469 | if (map_proj_id .eq. 3) goto 2000 |
---|
470 | if (map_proj_id .eq. 4) goto 2000 |
---|
471 | end if |
---|
472 | goto 5000 |
---|
473 | |
---|
474 | ! |
---|
475 | ! use map projection routines in map_proj_wrf_subaa.f |
---|
476 | ! |
---|
477 | 2000 continue |
---|
478 | |
---|
479 | call latlon_to_ij( y_lat, x_lon, dumxi, dumyj ) |
---|
480 | x_i = dumxi - 1.0 |
---|
481 | y_j = dumyj - 1.0 |
---|
482 | goto 8000 |
---|
483 | |
---|
484 | |
---|
485 | ! |
---|
486 | ! do it by search/minimization of distance, |
---|
487 | ! using the outer grid xlon2d/ylat2d values from WRF met file |
---|
488 | ! |
---|
489 | 5000 continue |
---|
490 | |
---|
491 | ! |
---|
492 | ! first locate the i,j for which the lon,lat at |
---|
493 | ! i+0.5,j+0.5 are closest to x_lon, y_lat |
---|
494 | ! |
---|
495 | dumr2min = 1.0e30 |
---|
496 | dumcos = cos( y_lat ) |
---|
497 | do j = 0, ny-1 |
---|
498 | do i = 0, nx-1 |
---|
499 | dumlat = 0.25*( ylat2d(i,j ) + ylat2d(i+1,j ) + & |
---|
500 | ylat2d(i,j+1) + ylat2d(i+1,j+1) ) |
---|
501 | dumlon = 0.25*( xlon2d(i,j ) + xlon2d(i+1,j ) + & |
---|
502 | xlon2d(i,j+1) + xlon2d(i+1,j+1) ) |
---|
503 | dumr2 = (y_lat-dumlat)**2 + ((x_lon-dumlon)*dumcos)**2 |
---|
504 | if (dumr2 .lt. dumr2min) then |
---|
505 | dumr2min = dumr2 |
---|
506 | ib = i |
---|
507 | jb = j |
---|
508 | end if |
---|
509 | end do |
---|
510 | end do |
---|
511 | i = ib |
---|
512 | j = jb |
---|
513 | ip = i+1 |
---|
514 | jp = j+1 |
---|
515 | |
---|
516 | ! |
---|
517 | ! next determine the position between i/i+1 & j/j+1 |
---|
518 | ! that is closest x_lon, y_lat |
---|
519 | ! |
---|
520 | xxyydelmin = 5.0e-8*max( abs(x_lon), abs(y_lat) ) |
---|
521 | xxyydel = 1.00 |
---|
522 | xxcen = 0.5 |
---|
523 | yycen = 0.5 |
---|
524 | iijjlohi = 3 |
---|
525 | ipass = 0 |
---|
526 | |
---|
527 | ! write(*,9510) |
---|
528 | ! write(*,9520) 0, (i+xxcen), (j+yycen), x_lon, y_lat |
---|
529 | 9510 format( / 'll_to_xyindex' ) |
---|
530 | 9520 format( 'ipass, x_i, y_j, lon, lat', i3, 4f14.7 ) |
---|
531 | |
---|
532 | 5200 ipass = ipass + 1 |
---|
533 | dumr2min = 1.0e30 |
---|
534 | do jj = -iijjlohi, iijjlohi |
---|
535 | do ii = -iijjlohi, iijjlohi |
---|
536 | fhi = xxcen + ii*xxyydel |
---|
537 | ghi = yycen + jj*xxyydel |
---|
538 | flo = 1.0 - fhi |
---|
539 | glo = 1.0 - ghi |
---|
540 | dumlat = glo*(flo*ylat2d(i,j ) + fhi*ylat2d(ip,j )) + & |
---|
541 | ghi*(flo*ylat2d(i,jp) + fhi*ylat2d(ip,jp)) |
---|
542 | dumlon = glo*(flo*xlon2d(i,j ) + fhi*xlon2d(ip,j )) + & |
---|
543 | ghi*(flo*xlon2d(i,jp) + fhi*xlon2d(ip,jp)) |
---|
544 | dumr2 = (y_lat-dumlat)**2 + ((x_lon-dumlon)*dumcos)**2 |
---|
545 | if (dumr2 .lt. dumr2min) then |
---|
546 | dumr2min = dumr2 |
---|
547 | xxcenb = fhi |
---|
548 | yycenb = ghi |
---|
549 | end if |
---|
550 | end do |
---|
551 | end do |
---|
552 | xxcen = xxcenb |
---|
553 | yycen = yycenb |
---|
554 | ! write(*,9520) ipass, (i+xxcen), (j+yycen), dumlon, dumlat |
---|
555 | if (xxyydel .gt. xxyydelmin) then |
---|
556 | xxyydel = xxyydel*0.5 |
---|
557 | if (ipass .eq. 4) iijjlohi = 2 |
---|
558 | goto 5200 |
---|
559 | end if |
---|
560 | ! write(*,9520) ipass, (i+xxcen), (j+yycen), dumlon, dumlat |
---|
561 | |
---|
562 | x_i = i + xxcen |
---|
563 | y_j = j + yycen |
---|
564 | |
---|
565 | |
---|
566 | ! |
---|
567 | ! check for x_i, y_j in bounds before returning |
---|
568 | ! |
---|
569 | 8000 continue |
---|
570 | if ((x_i .lt. -grace) .or. (x_i .gt. grace+real(nx-1)) .or. & |
---|
571 | (y_j .lt. -grace) .or. (y_j .gt. grace+real(ny-1))) then |
---|
572 | write(*,'(/a/a,1p,2e12.3/a,1p,2e12.3)') & |
---|
573 | '*** ll_to_xyindex_wrf error -- x_i, y_j out of bounds', & |
---|
574 | ' x_lon, y_lat =', x_lon_sv, y_lat_sv, & |
---|
575 | ' x_i, y_j =', x_i, y_j |
---|
576 | ! pause |
---|
577 | ! ylat2d(5000,2000)=0. |
---|
578 | ! stop |
---|
579 | end if |
---|
580 | return |
---|
581 | |
---|
582 | end subroutine ll_to_xyindex_wrf |
---|
583 | |
---|
584 | |
---|
585 | !----------------------------------------------------------------------- |
---|
586 | subroutine ll_to_xymeter_wrf( x_lon, y_lat, xmeter, ymeter ) |
---|
587 | ! |
---|
588 | ! calculates xy "grid meter" coordinates from longitude/latitude |
---|
589 | ! |
---|
590 | ! arguments |
---|
591 | ! x_lon, y_lat - input - longitude and latitude in degrees |
---|
592 | ! xmeter, ymeter - output = "grid meter" coordinates on the mother grid. |
---|
593 | ! |
---|
594 | ! include 'includepar' |
---|
595 | ! include 'includecom' |
---|
596 | use par_mod |
---|
597 | use com_mod |
---|
598 | implicit none |
---|
599 | |
---|
600 | ! arguments |
---|
601 | real :: xmeter, ymeter, x_lon, y_lat |
---|
602 | |
---|
603 | ! local variables |
---|
604 | real :: x_i, y_j |
---|
605 | |
---|
606 | call ll_to_xyindex_wrf( x_lon, y_lat, x_i, y_j ) |
---|
607 | xmeter = xmet0 + dx*x_i |
---|
608 | ymeter = ymet0 + dy*y_j |
---|
609 | |
---|
610 | return |
---|
611 | end subroutine ll_to_xymeter_wrf |
---|
612 | |
---|
613 | |
---|
614 | !----------------------------------------------------------------------- |
---|
615 | subroutine test_xyindex_to_ll_wrf( lgrid ) |
---|
616 | ! |
---|
617 | ! tests the map projection routines by comparing |
---|
618 | ! lat,lon from projection routine against the |
---|
619 | ! lat,lon from the WRF met. files |
---|
620 | ! |
---|
621 | ! arguments |
---|
622 | ! lgrid - input - grid identifier |
---|
623 | ! |
---|
624 | ! include 'includepar' |
---|
625 | ! include 'includecom' |
---|
626 | use par_mod |
---|
627 | use com_mod |
---|
628 | |
---|
629 | ! arguments |
---|
630 | implicit none |
---|
631 | integer :: lgrid |
---|
632 | ! local variables |
---|
633 | integer :: idum, ix, jdum, jy |
---|
634 | integer :: map_set_proj_code |
---|
635 | real :: dumdx, dumxi, dumyj |
---|
636 | real :: dumlat, dumlatb, dumlon, dumlonb |
---|
637 | real :: err |
---|
638 | real :: rmserr |
---|
639 | |
---|
640 | |
---|
641 | |
---|
642 | ! |
---|
643 | ! check if map projection is lambert conformal or polar stereographic |
---|
644 | ! |
---|
645 | if (map_proj_id .eq. 1) then |
---|
646 | map_set_proj_code = 3 ! lambert conformal |
---|
647 | else if (map_proj_id .eq. 2) then |
---|
648 | map_set_proj_code = 5 ! polar stereographic |
---|
649 | else if (map_proj_id .eq. 3) then |
---|
650 | map_set_proj_code = 1 ! mercator |
---|
651 | else if (map_proj_id .eq. 4) then |
---|
652 | map_set_proj_code = 0 ! lat/lon |
---|
653 | else |
---|
654 | write(*,'(/ 10(a/) )') & |
---|
655 | '************************************************************', & |
---|
656 | '* *', & |
---|
657 | '* WARNING - map projection is not polar sterographic *', & |
---|
658 | '* or lambert conformal *', & |
---|
659 | '* *', & |
---|
660 | '* x,y <--> lat,lon conversions will be done *', & |
---|
661 | '* by interpolation & searching, and will *', & |
---|
662 | '* have limited accuracy near poles *', & |
---|
663 | '* *', & |
---|
664 | '************************************************************' |
---|
665 | map_proj_method = 0 |
---|
666 | return |
---|
667 | end if |
---|
668 | if (lgrid .gt. 0) goto 2000 |
---|
669 | |
---|
670 | ! |
---|
671 | ! make call to map projection setup routine |
---|
672 | ! |
---|
673 | ! (The 0.999812 factor is due to different earth_radius |
---|
674 | ! values in wrfsi code (6370.0 vs 6371.2).) |
---|
675 | dumdx = dx*0.999812 ! WA 2/11/10 is this right? |
---|
676 | ! dumdx = dx*0.9987242 ! WA 2/11/10 is this right? |
---|
677 | ! dumdx = dx*0.9979203 ! WA 2/11/10 is this right? |
---|
678 | coefdx=0.999812 |
---|
679 | ! write(*,'(/2a,2i5)') & |
---|
680 | ! 'test_xyindex_to_ll_wrf calling map_set -- ', & |
---|
681 | ! 'map_proj_id, map_set_proj_code =', & |
---|
682 | ! map_proj_id, map_set_proj_code |
---|
683 | call map_set( map_set_proj_code, & |
---|
684 | ylat2d(0,0), xlon2d(0,0), dumdx, & |
---|
685 | map_stdlon, map_truelat1, map_truelat2, & |
---|
686 | nx, ny ) |
---|
687 | map_proj_method = 1 |
---|
688 | |
---|
689 | ! |
---|
690 | ! compute lat,lon from xi,jy at 9 points (center, 4 corners, |
---|
691 | ! and 4 boundary midpoints) |
---|
692 | ! compare to lon,lat read from the WRF met. file and report rmserr |
---|
693 | ! |
---|
694 | 2000 continue |
---|
695 | rmserr = 0.0 |
---|
696 | do jdum = 0, 2 |
---|
697 | do idum = 0, 2 |
---|
698 | if (lgrid .le. 0) then |
---|
699 | jy = nint( 0.5*real((ny-1)*jdum) ) |
---|
700 | ix = nint( 0.5*real((nx-1)*idum) ) |
---|
701 | dumyj = 1 + jy |
---|
702 | dumxi = 1 + ix |
---|
703 | dumlatb = ylat2d(ix,jy) |
---|
704 | dumlonb = xlon2d(ix,jy) |
---|
705 | else |
---|
706 | jy = nint( 0.5*real((nyn(lgrid)-1)*jdum) ) |
---|
707 | ix = nint( 0.5*real((nxn(lgrid)-1)*idum) ) |
---|
708 | dumyj = 1.0 + (jy/yresoln(lgrid)) + yln(lgrid) |
---|
709 | dumxi = 1.0 + (ix/xresoln(lgrid)) + xln(lgrid) |
---|
710 | dumlatb = ylat2dn(ix,jy,lgrid) |
---|
711 | dumlonb = xlon2dn(ix,jy,lgrid) |
---|
712 | end if |
---|
713 | call ij_to_latlon( dumxi, dumyj, dumlat, dumlon ) |
---|
714 | err = (dumlat-dumlatb)**2 + & |
---|
715 | ((dumlon-dumlonb)*cos(dumlatb))**2 |
---|
716 | ! print*,'err',idum,jdum,dumlat,dumlatb,dumlon,dumlonb |
---|
717 | rmserr = rmserr + err |
---|
718 | ! print*,dumxi, dumyj, dumlat, dumlon,dumlatb,dumlonb |
---|
719 | end do |
---|
720 | end do |
---|
721 | rmserr = 111.2*sqrt( rmserr/9.0 ) |
---|
722 | ! print*,'rmserr',rmserr |
---|
723 | if (rmserr .le. 0.02*(1.0e-3*dumdx)) then |
---|
724 | write(*,'(/a,i3,1pe10.2)') & |
---|
725 | 'test_xyindex_to_ll_wrf -- lgrid, rmserr (km) =', lgrid, rmserr |
---|
726 | endif |
---|
727 | ! |
---|
728 | ! if rms error exceeds 0.02*dx, something is wrong |
---|
729 | ! do not use the map projection routines in this case |
---|
730 | ! |
---|
731 | if (rmserr .gt. 0.02*(1.0e-3*dumdx)) then |
---|
732 | ! print*,'pb with rmserr. try an other earth radius' |
---|
733 | dumdx = dx ! try another earth radius |
---|
734 | coefdx=1. |
---|
735 | |
---|
736 | write(*,'(/2a,2i5)') & |
---|
737 | 'test_xyindex_to_ll_wrf calling map_set -- ', & |
---|
738 | 'map_proj_id, map_set_proj_code =', & |
---|
739 | map_proj_id, map_set_proj_code |
---|
740 | call map_set( map_set_proj_code, & |
---|
741 | ylat2d(0,0), xlon2d(0,0), dumdx, & |
---|
742 | map_stdlon, map_truelat1, map_truelat2, & |
---|
743 | nx, ny ) |
---|
744 | map_proj_method = 1 |
---|
745 | |
---|
746 | ! |
---|
747 | ! compute lat,lon from xi,jy at 9 points (center, 4 corners, |
---|
748 | ! and 4 boundary midpoints) |
---|
749 | ! compare to lon,lat read from the WRF met. file and report rmserr |
---|
750 | ! |
---|
751 | rmserr = 0.0 |
---|
752 | do jdum = 0, 2 |
---|
753 | do idum = 0, 2 |
---|
754 | if (lgrid .le. 0) then |
---|
755 | jy = nint( 0.5*real((ny-1)*jdum) ) |
---|
756 | ix = nint( 0.5*real((nx-1)*idum) ) |
---|
757 | dumyj = 1 + jy |
---|
758 | dumxi = 1 + ix |
---|
759 | dumlatb = ylat2d(ix,jy) |
---|
760 | dumlonb = xlon2d(ix,jy) |
---|
761 | else |
---|
762 | jy = nint( 0.5*real((nyn(lgrid)-1)*jdum) ) |
---|
763 | ix = nint( 0.5*real((nxn(lgrid)-1)*idum) ) |
---|
764 | dumyj = 1.0 + (jy/yresoln(lgrid)) + yln(lgrid) |
---|
765 | dumxi = 1.0 + (ix/xresoln(lgrid)) + xln(lgrid) |
---|
766 | dumlatb = ylat2dn(ix,jy,lgrid) |
---|
767 | dumlonb = xlon2dn(ix,jy,lgrid) |
---|
768 | end if |
---|
769 | call ij_to_latlon( dumxi, dumyj, dumlat, dumlon ) |
---|
770 | err = (dumlat-dumlatb)**2 + & |
---|
771 | ((dumlon-dumlonb)*cos(dumlatb))**2 |
---|
772 | rmserr = rmserr + err |
---|
773 | end do |
---|
774 | end do |
---|
775 | rmserr = 111.2*sqrt( rmserr/9.0 ) |
---|
776 | |
---|
777 | write(*,'(/a,i3,1pe10.2)') & |
---|
778 | 'test_xyindex_to_ll_wrf -- lgrid, rmserr (km) =', lgrid, rmserr |
---|
779 | |
---|
780 | |
---|
781 | endif |
---|
782 | |
---|
783 | ! |
---|
784 | ! if rms error exceeds 0.02*dx, something is wrong |
---|
785 | ! do not use the map projection routines in this case |
---|
786 | ! |
---|
787 | dumdx = dx |
---|
788 | if (lgrid .gt. 0) dumdx = dxn(lgrid) |
---|
789 | if (rmserr .gt. 0.02*(1.0e-3*dumdx)) then |
---|
790 | map_proj_method = 0 |
---|
791 | write(*,'(/ 9(a/) )') & |
---|
792 | '************************************************************', & |
---|
793 | '* *', & |
---|
794 | '* WARNING - the coordinate transfo error exceeds 0.02*dx *', & |
---|
795 | '* *', & |
---|
796 | '************************************************************' |
---|
797 | ! print*,'projection set to',map_proj_method |
---|
798 | print*,'problem with the projection. wrf+flexpart stops' |
---|
799 | stop |
---|
800 | end if |
---|
801 | |
---|
802 | return |
---|
803 | end subroutine test_xyindex_to_ll_wrf |
---|
804 | |
---|
805 | |
---|