1 | ! SUPPORTED PROJECTIONS |
---|
2 | ! --------------------- |
---|
3 | ! Cylindrical Lat/Lon (code = PROJ_LATLON) |
---|
4 | ! Mercator (code = PROJ_MERC) |
---|
5 | ! Lambert Conformal (code = PROJ_LC) |
---|
6 | ! Polar Stereographic (code = PROJ_PS) |
---|
7 | ! |
---|
8 | ! REMARKS |
---|
9 | ! ------- |
---|
10 | ! The routines contained within were adapted from routines |
---|
11 | ! obtained from the NCEP w3 library. The original NCEP routines were |
---|
12 | ! less |
---|
13 | ! flexible (e.g., polar-stereo routines only supported truelat of |
---|
14 | ! 60N/60S) |
---|
15 | ! than what we needed, so modifications based on equations in Hoke, |
---|
16 | ! Hayes, and |
---|
17 | ! Renninger (AFGWC/TN/79-003) were added to improve the flexibility. |
---|
18 | ! Additionally, coding was improved to F90 standards and the routines |
---|
19 | ! were |
---|
20 | ! combined into this module. |
---|
21 | !----------------------------------------------------------------------- |
---|
22 | |
---|
23 | |
---|
24 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
25 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
26 | !dis |
---|
27 | !dis open source license/disclaimer, forecast systems laboratory |
---|
28 | !dis noaa/oar/fsl, 325 broadway boulder, co 80305 |
---|
29 | !dis |
---|
30 | !dis this software is distributed under the open source definition, |
---|
31 | !dis which may be found at http://www.opensource.org/osd.html. |
---|
32 | !dis |
---|
33 | !dis in particular, redistribution and use in source and binary forms, |
---|
34 | !dis with or without modification, are permitted provided that the |
---|
35 | !dis following conditions are met: |
---|
36 | !dis |
---|
37 | !dis - redistributions of source code must retain this notice, this |
---|
38 | !dis list of conditions and the following disclaimer. |
---|
39 | !dis |
---|
40 | !dis - redistributions in binary form must provide access to this |
---|
41 | !dis notice, this list of conditions and the following disclaimer, and |
---|
42 | !dis the underlying source code. |
---|
43 | !dis |
---|
44 | !dis - all modifications to this software must be clearly documented, |
---|
45 | !dis and are solely the responsibility of the agent making the |
---|
46 | !dis modifications. |
---|
47 | !dis |
---|
48 | !dis - if significant modifications or enhancements are made to this |
---|
49 | !dis software, the fsl software policy manager |
---|
50 | !dis (softwaremgr@fsl.noaa.gov) should be notified. |
---|
51 | !dis |
---|
52 | !dis this software and its documentation are in the public domain |
---|
53 | !dis and are furnished "as is." the authors, the united states |
---|
54 | !dis government, its instrumentalities, officers, employees, and |
---|
55 | !dis agents make no warranty, express or implied, as to the usefulness |
---|
56 | !dis of the software and documentation for any purpose. they assume |
---|
57 | !dis no responsibility (1) for the use of the software and |
---|
58 | !dis documentation; or (2) to provide technical support to users. |
---|
59 | !dis |
---|
60 | !dis |
---|
61 | |
---|
62 | ! module that defines constants, data structures, and |
---|
63 | ! routines used to convert grid indices to lat/lon |
---|
64 | ! and vice versa. |
---|
65 | ! |
---|
66 | ! supported projections |
---|
67 | ! --------------------- |
---|
68 | ! cylindrical lat/lon (code = proj_latlon) |
---|
69 | ! mercator (code = proj_merc) |
---|
70 | ! lambert conformal (code = proj_lc) |
---|
71 | ! polar stereographic (code = proj_ps) |
---|
72 | ! |
---|
73 | ! remarks |
---|
74 | ! ------- |
---|
75 | ! the routines contained within were adapted from routines |
---|
76 | ! obtained from the ncep w3 library. the original ncep routines were less |
---|
77 | ! flexible (e.g., polar-stereo routines only supported truelat of 60n/60s) |
---|
78 | ! than what we needed, so modifications based on equations in hoke, hayes, and |
---|
79 | ! renninger (afgwc/tn/79-003) were added to improve the flexibility. |
---|
80 | ! additionally, coding was improved to f90 standards and the routines were |
---|
81 | ! combined into this module. |
---|
82 | ! |
---|
83 | ! assumptions |
---|
84 | ! ----------- |
---|
85 | ! grid definition: |
---|
86 | ! for mercator, lambert conformal, and polar-stereographic projections, |
---|
87 | ! the routines within assume the following: |
---|
88 | ! |
---|
89 | ! 1. grid is dimensioned (i,j) where i is the east-west direction, |
---|
90 | ! positive toward the east, and j is the north-south direction, |
---|
91 | ! positive toward the north. |
---|
92 | ! 2. origin is at (1,1) and is located at the southwest corner, |
---|
93 | ! regardless of hemispere. |
---|
94 | ! 3. grid spacing (dx) is always positive. |
---|
95 | ! 4. values of true latitudes must be positive for nh domains |
---|
96 | ! and negative for sh domains. |
---|
97 | ! |
---|
98 | ! for the latlon projection, the grid origin may be at any of the |
---|
99 | ! corners, and the deltalat and deltalon values can be signed to |
---|
100 | ! account for this using the following convention: |
---|
101 | ! origin location deltalat sign deltalon sign |
---|
102 | ! --------------- ------------- ------------- |
---|
103 | ! sw corner + + |
---|
104 | ! ne corner - - |
---|
105 | ! nw corner - + |
---|
106 | ! se corner + - |
---|
107 | ! |
---|
108 | ! data definitions: |
---|
109 | ! 1. any arguments that are a latitude value are expressed in |
---|
110 | ! degrees north with a valid range of -90 -> 90 |
---|
111 | ! 2. any arguments that are a longitude value are expressed in |
---|
112 | ! degrees east with a valid range of -180 -> 180. |
---|
113 | ! 3. distances are in meters and are always positive. |
---|
114 | ! 4. the standard longitude (stdlon) is defined as the longitude |
---|
115 | ! line which is parallel to the y-axis (j-direction), along |
---|
116 | ! which latitude increases (not the absolute value of latitude, but |
---|
117 | ! the actual latitude, such that latitude increases continuously |
---|
118 | ! from the south pole to the north pole) as j increases. |
---|
119 | ! 5. one true latitude value is required for polar-stereographic and |
---|
120 | ! mercator projections, and defines at which latitude the |
---|
121 | ! grid spacing is true. for lambert conformal, two true latitude |
---|
122 | ! values must be specified, but may be set equal to each other to |
---|
123 | ! specify a tangent projection instead of a secant projection. |
---|
124 | ! |
---|
125 | ! usage |
---|
126 | ! ----- |
---|
127 | ! to use the routines in this module, the calling routines must have the |
---|
128 | ! following statement at the beginning of its declaration block: |
---|
129 | ! use map_utils |
---|
130 | ! |
---|
131 | ! the use of the module not only provides access to the necessary routines, |
---|
132 | ! but also defines a structure of type (proj_info) that can be used |
---|
133 | ! to declare a variable of the same type to hold your map projection |
---|
134 | ! information. it also defines some integer parameters that contain |
---|
135 | ! the projection codes so one only has to use those variable names rather |
---|
136 | ! than remembering the acutal code when using them. the basic steps are |
---|
137 | ! as follows: |
---|
138 | ! |
---|
139 | ! 1. ensure the "use map_utils" is in your declarations. |
---|
140 | ! 2. declare the projection information structure as type(proj_info): |
---|
141 | ! type(proj_info) :: proj |
---|
142 | ! 3. populate your structure by calling the map_set routine: |
---|
143 | ! call map_set(code,lat1,lon1,dx,stdlon,truelat1,truelat2,nx,ny,proj) |
---|
144 | ! where: |
---|
145 | ! code (input) = one of proj_latlon, proj_merc, proj_lc, or proj_ps |
---|
146 | ! lat1 (input) = latitude of grid origin point (i,j)=(1,1) |
---|
147 | ! (see assumptions!) |
---|
148 | ! lon1 (input) = longitude of grid origin |
---|
149 | ! dx (input) = grid spacing in meters (ignored for latlon projections) |
---|
150 | ! stdlon (input) = standard longitude for proj_ps and proj_lc, |
---|
151 | ! deltalon (see assumptions) for proj_latlon, |
---|
152 | ! ignored for proj_merc |
---|
153 | ! truelat1 (input) = 1st true latitude for proj_ps, proj_lc, and |
---|
154 | ! proj_merc, deltalat (see assumptions) for proj_latlon |
---|
155 | ! truelat2 (input) = 2nd true latitude for proj_lc, |
---|
156 | ! ignored for all others. |
---|
157 | ! nx = number of points in east-west direction |
---|
158 | ! ny = number of points in north-south direction |
---|
159 | ! proj (output) = the structure of type (proj_info) that will be fully |
---|
160 | ! populated after this call |
---|
161 | ! |
---|
162 | ! 4. now that the proj structure is populated, you may call any |
---|
163 | ! of the following routines: |
---|
164 | ! |
---|
165 | ! latlon_to_ij(proj, lat, lon, i, j) |
---|
166 | ! ij_to_latlon(proj, i, j, lat, lon) |
---|
167 | ! truewind_to_gridwind(lon, proj, ugrid, vgrid, utrue, vtrue) |
---|
168 | ! gridwind_to_truewind(lon, proj, utrue, vtrue, ugrid, vgrid) |
---|
169 | ! compare_projections(proj1, proj2, same_proj) |
---|
170 | ! |
---|
171 | ! it is incumbent upon the calling routine to determine whether or |
---|
172 | ! not the values returned are within your domain bounds. all values |
---|
173 | ! of i, j, lat, and lon are real values. |
---|
174 | ! |
---|
175 | ! |
---|
176 | ! references |
---|
177 | ! ---------- |
---|
178 | ! hoke, hayes, and renninger, "map preojections and grid systems for |
---|
179 | ! meteorological applications." afgwc/tn-79/003(rev), air weather |
---|
180 | ! service, 1985. |
---|
181 | ! |
---|
182 | ! ncar mm5v3 modeling system, regridder program, module_first_guess_map.f |
---|
183 | ! ncep routines w3fb06, w3fb07, w3fb08, w3fb09, w3fb11, w3fb12 |
---|
184 | ! |
---|
185 | ! history |
---|
186 | ! ------- |
---|
187 | ! 27 mar 2001 - original version |
---|
188 | ! brent l. shaw, noaa/fsl (csu/cira) |
---|
189 | ! 02 apr 2001 - added routines to rotate winds from true to grid |
---|
190 | ! and vice versa. |
---|
191 | ! brent l. shaw, noaa/fsl (csu/cira) |
---|
192 | ! 09 apr 2001 - added compare_projections routine to compare two |
---|
193 | ! sets of projection parameters. |
---|
194 | ! |
---|
195 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
196 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
197 | |
---|
198 | |
---|
199 | |
---|
200 | |
---|
201 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
202 | ! ! define data structures to define various projections |
---|
203 | ! |
---|
204 | ! type proj_info |
---|
205 | ! |
---|
206 | ! logical proj_init ! flag to indicate if this struct is |
---|
207 | ! ! ready for use |
---|
208 | ! logical proj_cyclic ! flag indicating if this grid |
---|
209 | ! ! is cyclic in the longitudinal |
---|
210 | ! ! direction...happens with |
---|
211 | ! ! global lat/lon grids like gfs/avn |
---|
212 | ! integer proj_code ! integer code for projection type |
---|
213 | ! integer proj_nx |
---|
214 | ! integer proj_ny |
---|
215 | ! real proj_lat1 ! sw latitude (1,1) in degrees (-90->90n) |
---|
216 | ! real proj_lon1 ! sw longitude (1,1) in degrees (-180->180e) |
---|
217 | ! real proj_dx ! grid spacing in meters at truelats, used |
---|
218 | ! ! only for ps, lc, and merc projections |
---|
219 | ! real proj_dlat ! lat increment for lat/lon grids |
---|
220 | ! real proj_dlon ! lon increment for lat/lon grids |
---|
221 | ! real proj_clat ! center latitude of grid |
---|
222 | ! real proj_clon ! center longitude of grid |
---|
223 | ! real proj_stdlon ! longitude parallel to y-axis (-180->180e) |
---|
224 | ! real proj_truelat1 ! first true latitude (all projections) |
---|
225 | ! real proj_truelat2 ! second true lat (lc only) |
---|
226 | ! real proj_hemi ! 1 for nh, -1 for sh |
---|
227 | ! real proj_cone ! cone factor for lc projections |
---|
228 | ! real proj_polei ! computed i-location of pole point |
---|
229 | ! real proj_polej ! computed j-location of pole point |
---|
230 | ! real proj_rsw ! computed radius to sw corner |
---|
231 | ! real proj_rebydx ! earth radius divided by dx |
---|
232 | ! |
---|
233 | ! end type proj_info |
---|
234 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
235 | |
---|
236 | |
---|
237 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
238 | subroutine map_init |
---|
239 | ! initializes the map projection structure to missing values |
---|
240 | |
---|
241 | use wrf_map_utils_mod |
---|
242 | implicit none |
---|
243 | |
---|
244 | ! include 'include_wrf_map_utils' |
---|
245 | pi = 3.1415927 |
---|
246 | deg_per_rad = 180./pi |
---|
247 | rad_per_deg = pi / 180. |
---|
248 | |
---|
249 | ! mean earth radius in m. the value below is consistent |
---|
250 | ! with nceps routines and grids. |
---|
251 | earth_radius_m = 6370000. |
---|
252 | ! earth_radius_m = 6371200. |
---|
253 | ! earth_radius_m = 6368750. |
---|
254 | |
---|
255 | proj_latlon = 0 |
---|
256 | proj_merc = 1 |
---|
257 | proj_lc = 3 |
---|
258 | proj_ps = 5 |
---|
259 | proj_rotlat = 203 |
---|
260 | |
---|
261 | proj_lat1 = -999.9 |
---|
262 | proj_lon1 = -999.9 |
---|
263 | proj_dx = -999.9 |
---|
264 | proj_stdlon = -999.9 |
---|
265 | proj_truelat1 = -999.9 |
---|
266 | proj_truelat2 = -999.9 |
---|
267 | proj_hemi = 0.0 |
---|
268 | proj_cone = -999.9 |
---|
269 | proj_polei = -999.9 |
---|
270 | proj_polej = -999.9 |
---|
271 | proj_rsw = -999.9 |
---|
272 | proj_init = .false. |
---|
273 | proj_nx = -99 |
---|
274 | proj_ny = -99 |
---|
275 | proj_cyclic = .false. |
---|
276 | |
---|
277 | return |
---|
278 | end subroutine map_init |
---|
279 | |
---|
280 | |
---|
281 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
282 | subroutine map_set( & |
---|
283 | proj_code_in, lat1, lon1, dx, & |
---|
284 | stdlon, truelat1, truelat2, & |
---|
285 | idim, jdim ) |
---|
286 | ! given a partially filled proj_info structure, this routine computes |
---|
287 | ! polei, polej, rsw, and cone (if lc projection) to complete the |
---|
288 | ! structure. this allows us to eliminate redundant calculations when |
---|
289 | ! calling the coordinate conversion routines multiple times for the |
---|
290 | ! same map. |
---|
291 | ! this will generally be the first routine called when a user wants |
---|
292 | ! to be able to use the coordinate conversion routines, and it |
---|
293 | ! will call the appropriate routines based on the |
---|
294 | ! proj_code which indicates which projection type this is. |
---|
295 | use wrf_map_utils_mod |
---|
296 | |
---|
297 | implicit none |
---|
298 | |
---|
299 | ! declare arguments |
---|
300 | integer :: proj_code_in |
---|
301 | real :: lat1 |
---|
302 | real :: lon1 |
---|
303 | real :: dx |
---|
304 | real :: stdlon |
---|
305 | real :: truelat1 |
---|
306 | real :: truelat2 |
---|
307 | integer :: idim |
---|
308 | integer :: jdim |
---|
309 | |
---|
310 | ! local variables |
---|
311 | real :: center_i,center_j |
---|
312 | real :: center_lat, center_lon |
---|
313 | |
---|
314 | ! include 'include_wrf_map_utils' |
---|
315 | |
---|
316 | ! executable code |
---|
317 | |
---|
318 | proj_code = proj_code_in |
---|
319 | |
---|
320 | ! first, check for validity of mandatory variables in proj |
---|
321 | if ( abs(lat1) .gt. 90.001 ) then |
---|
322 | print '(a)', 'latitude of origin corner required as follows:' |
---|
323 | print '(a)', ' -90n <= lat1 < = 90.n' |
---|
324 | stop 'map_init' |
---|
325 | endif |
---|
326 | if ( abs(lon1) .gt. 180.) then |
---|
327 | print '(a)', 'longitude of origin required as follows:' |
---|
328 | print '(a)', ' -180e <= lon1 <= 180w' |
---|
329 | stop 'map_init' |
---|
330 | endif |
---|
331 | if ((dx .le. 0.).and.(proj_code .ne. proj_latlon)) then |
---|
332 | print '(a)', 'require grid spacing (dx) in meters be positive!' |
---|
333 | stop 'map_init' |
---|
334 | endif |
---|
335 | if ((abs(stdlon) .gt. 180.).and.(proj_code .ne. proj_merc)) then |
---|
336 | print '(a)', 'need orientation longitude (stdlon) as: ' |
---|
337 | print '(a)', ' -180e <= lon1 <= 180w' |
---|
338 | stop 'map_init' |
---|
339 | endif |
---|
340 | if (abs(truelat1).gt.90.) then |
---|
341 | print '(a)', 'set true latitude 1 for all projections!' |
---|
342 | stop 'map_init' |
---|
343 | endif |
---|
344 | |
---|
345 | call map_init |
---|
346 | proj_code = proj_code_in |
---|
347 | proj_lat1 = lat1 |
---|
348 | proj_lon1 = lon1 |
---|
349 | proj_dx = dx |
---|
350 | proj_stdlon = stdlon |
---|
351 | proj_truelat1 = truelat1 |
---|
352 | proj_truelat2 = truelat2 |
---|
353 | proj_nx = idim |
---|
354 | proj_ny = jdim |
---|
355 | if (proj_code .ne. proj_latlon) then |
---|
356 | proj_dx = dx |
---|
357 | if (truelat1 .lt. 0.) then |
---|
358 | proj_hemi = -1.0 |
---|
359 | else |
---|
360 | proj_hemi = 1.0 |
---|
361 | endif |
---|
362 | proj_rebydx = earth_radius_m / dx |
---|
363 | endif |
---|
364 | |
---|
365 | |
---|
366 | if (proj_code .eq. proj_ps) then |
---|
367 | !print '(a)', 'setting up polar stereographic map...' |
---|
368 | call set_ps |
---|
369 | |
---|
370 | else if (proj_code .eq. proj_lc) then |
---|
371 | !print '(a)', 'setting up lambert conformal map...' |
---|
372 | if (abs(proj_truelat2) .gt. 90.) then |
---|
373 | print '(a)', & |
---|
374 | 'second true latitude not set, assuming a tangent' |
---|
375 | print '(a,f10.3)', 'projection at truelat1: ', proj_truelat1 |
---|
376 | proj_truelat2=proj_truelat1 |
---|
377 | else |
---|
378 | ! ensure truelat1 < truelat2 |
---|
379 | proj_truelat1 = min(truelat1,truelat2) |
---|
380 | proj_truelat2 = max(truelat1,truelat2) |
---|
381 | endif |
---|
382 | call set_lc |
---|
383 | |
---|
384 | else if (proj_code .eq. proj_merc) then |
---|
385 | !print '(a)', 'setting up mercator map...' |
---|
386 | call set_merc |
---|
387 | |
---|
388 | else if (proj_code .eq. proj_latlon) then |
---|
389 | !print '(a)', 'setting up cylindrical equidistant latlon map...' |
---|
390 | ! convert lon1 to 0->360 notation |
---|
391 | if (proj_lon1 .lt. 0.) proj_lon1 = proj_lon1 + 360. |
---|
392 | proj_dlat = truelat1 |
---|
393 | proj_dlon = stdlon |
---|
394 | if (nint(proj_dlon*real(proj_nx)) .eq. 360) proj_cyclic = .true. |
---|
395 | |
---|
396 | else |
---|
397 | print '(a,i2)', 'unknown projection code: ', proj_code |
---|
398 | stop 'map_init' |
---|
399 | |
---|
400 | end if |
---|
401 | |
---|
402 | proj_init = .true. |
---|
403 | |
---|
404 | center_i = real(proj_nx+1)*0.5 |
---|
405 | center_j = real(proj_ny+1)*0.5 |
---|
406 | call ij_to_latlon(center_i,center_j,center_lat,center_lon) |
---|
407 | proj_clat = center_lat |
---|
408 | proj_clon = center_lon |
---|
409 | |
---|
410 | return |
---|
411 | end subroutine map_set |
---|
412 | |
---|
413 | |
---|
414 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
415 | subroutine latlon_to_ij(lat, lon, i, j) |
---|
416 | ! converts input lat/lon values to the cartesian (i,j) value |
---|
417 | ! for the given projection. |
---|
418 | |
---|
419 | use wrf_map_utils_mod |
---|
420 | implicit none |
---|
421 | |
---|
422 | ! arguments |
---|
423 | real :: lat |
---|
424 | real :: lon |
---|
425 | real :: i |
---|
426 | real :: j |
---|
427 | |
---|
428 | ! include 'include_wrf_map_utils' |
---|
429 | |
---|
430 | if (.not.proj_init) then |
---|
431 | print '(a)', 'you have not called map_set for this projection!' |
---|
432 | stop 'latlon_to_ij' |
---|
433 | endif |
---|
434 | |
---|
435 | if (proj_code .eq. proj_ps) then |
---|
436 | call llij_ps(lat,lon,i,j) |
---|
437 | |
---|
438 | else if (proj_code .eq. proj_lc) then |
---|
439 | call llij_lc(lat,lon,i,j) |
---|
440 | |
---|
441 | else if (proj_code .eq. proj_latlon) then |
---|
442 | call llij_latlon(lat,lon,i,j) |
---|
443 | ! |
---|
444 | else if (proj_code .eq. proj_merc) then |
---|
445 | call llij_merc(lat,lon,i,j) |
---|
446 | ! |
---|
447 | ! else if (proj_code .eq. proj_rotlat) then |
---|
448 | ! write(6,*) 'doing nothing in latlon_to_ij' |
---|
449 | |
---|
450 | else |
---|
451 | print '(a,i2)','unrecognized map projection code: ', proj_code |
---|
452 | stop 'latlon_to_ij' |
---|
453 | |
---|
454 | end if |
---|
455 | |
---|
456 | return |
---|
457 | end subroutine latlon_to_ij |
---|
458 | |
---|
459 | |
---|
460 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
461 | subroutine ij_to_latlon(i, j, lat, lon) |
---|
462 | ! computes geographical latitude and longitude for a given (i,j) point |
---|
463 | ! in a grid with a projection of proj |
---|
464 | |
---|
465 | use wrf_map_utils_mod |
---|
466 | implicit none |
---|
467 | |
---|
468 | ! arguments |
---|
469 | real :: i |
---|
470 | real :: j |
---|
471 | real :: lat |
---|
472 | real :: lon |
---|
473 | |
---|
474 | ! include 'include_wrf_map_utils' |
---|
475 | |
---|
476 | if (.not.proj_init) then |
---|
477 | print '(a)', 'you have not called map_set for this projection!' |
---|
478 | stop 'ij_to_latlon' |
---|
479 | endif |
---|
480 | |
---|
481 | if (proj_code .eq. proj_ps) then |
---|
482 | call ijll_ps(i, j, lat, lon) |
---|
483 | |
---|
484 | else if (proj_code .eq. proj_lc) then |
---|
485 | call ijll_lc(i, j, lat, lon) |
---|
486 | |
---|
487 | else if (proj_code .eq. proj_latlon) then |
---|
488 | call ijll_latlon(i, j, lat, lon) |
---|
489 | ! |
---|
490 | else if (proj_code .eq. proj_merc) then |
---|
491 | call ijll_merc(i, j, lat, lon) |
---|
492 | ! |
---|
493 | ! else if (proj_code .eq. proj_rotlat) then |
---|
494 | ! write(6,*) 'doing nothing in ij_to_latlon' |
---|
495 | |
---|
496 | else |
---|
497 | print '(a,i2)','unrecognized map projection code: ', proj_code |
---|
498 | stop 'ij_to_latlon' |
---|
499 | |
---|
500 | end if |
---|
501 | |
---|
502 | return |
---|
503 | end subroutine ij_to_latlon |
---|
504 | |
---|
505 | |
---|
506 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
507 | subroutine set_ps |
---|
508 | ! initializes a polar-stereographic map projection from the partially |
---|
509 | ! filled proj structure. this routine computes the radius to the |
---|
510 | ! southwest corner and computes the i/j location of the pole for use |
---|
511 | ! in llij_ps and ijll_ps. |
---|
512 | |
---|
513 | use wrf_map_utils_mod |
---|
514 | implicit none |
---|
515 | |
---|
516 | ! local vars |
---|
517 | real :: ala1 |
---|
518 | real :: alo1 |
---|
519 | real :: reflon |
---|
520 | real :: scale_top |
---|
521 | real :: dumlat, dumlon |
---|
522 | |
---|
523 | ! include 'include_wrf_map_utils' |
---|
524 | |
---|
525 | ! executable code |
---|
526 | reflon = proj_stdlon + 90. |
---|
527 | |
---|
528 | ! cone factor |
---|
529 | proj_cone = 1.0 |
---|
530 | |
---|
531 | ! compute numerator term of map scale factor |
---|
532 | scale_top = 1. + proj_hemi * sin(proj_truelat1 * rad_per_deg) |
---|
533 | |
---|
534 | ! compute radius to lower-left (sw) corner |
---|
535 | ala1 = proj_lat1 * rad_per_deg |
---|
536 | proj_rsw =proj_rebydx*cos(ala1)*scale_top/(1.+proj_hemi*sin(ala1)) |
---|
537 | |
---|
538 | ! find the pole point |
---|
539 | alo1 = (proj_lon1 - reflon) * rad_per_deg |
---|
540 | proj_polei = 1. - proj_rsw * cos(alo1) |
---|
541 | proj_polej = 1. - proj_hemi * proj_rsw * sin(alo1) |
---|
542 | ! print '(a,2f14.5)', 'set_ps - computed (i,j) of pole point: ', |
---|
543 | ! & proj_polei,proj_polej |
---|
544 | |
---|
545 | return |
---|
546 | end subroutine set_ps |
---|
547 | |
---|
548 | |
---|
549 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
550 | subroutine llij_ps( lat, lon, i, j ) |
---|
551 | ! given latitude (-90 to 90), longitude (-180 to 180), and the |
---|
552 | ! standard polar-stereographic projection information via the |
---|
553 | ! public proj structure, this routine returns the i/j indices which |
---|
554 | ! if within the domain range from 1->nx and 1->ny, respectively. |
---|
555 | use wrf_map_utils_mod |
---|
556 | |
---|
557 | implicit none |
---|
558 | |
---|
559 | ! delcare input arguments |
---|
560 | real :: lat |
---|
561 | real :: lon |
---|
562 | |
---|
563 | ! declare output arguments |
---|
564 | real :: i !(x-index) |
---|
565 | real :: j !(y-index) |
---|
566 | |
---|
567 | ! declare local variables |
---|
568 | real :: reflon |
---|
569 | real :: scale_top |
---|
570 | real :: ala |
---|
571 | real :: alo |
---|
572 | real :: rm |
---|
573 | |
---|
574 | ! include 'include_wrf_map_utils' |
---|
575 | |
---|
576 | ! begin code |
---|
577 | |
---|
578 | reflon = proj_stdlon + 90. |
---|
579 | |
---|
580 | ! compute numerator term of map scale factor |
---|
581 | |
---|
582 | scale_top = 1. + proj_hemi * sin(proj_truelat1 * rad_per_deg) |
---|
583 | |
---|
584 | ! find radius to desired point |
---|
585 | ala = lat * rad_per_deg |
---|
586 | rm = proj_rebydx * cos(ala) * scale_top/(1. + proj_hemi *sin(ala)) |
---|
587 | alo = (lon - reflon) * rad_per_deg |
---|
588 | i = proj_polei + rm * cos(alo) |
---|
589 | j = proj_polej + proj_hemi * rm * sin(alo) |
---|
590 | |
---|
591 | return |
---|
592 | end subroutine llij_ps |
---|
593 | |
---|
594 | |
---|
595 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
596 | subroutine ijll_ps( i, j, lat, lon ) |
---|
597 | |
---|
598 | ! this is the inverse routine of llij_ps. it returns the |
---|
599 | ! latitude and longitude of an i/j point given the projection info |
---|
600 | ! structure. |
---|
601 | use wrf_map_utils_mod |
---|
602 | |
---|
603 | implicit none |
---|
604 | |
---|
605 | ! declare input arguments |
---|
606 | real :: i ! column |
---|
607 | real :: j ! row |
---|
608 | |
---|
609 | ! declare output arguments |
---|
610 | real :: lat ! -90 -> 90 north |
---|
611 | real :: lon ! -180 -> 180 east |
---|
612 | |
---|
613 | ! local variables |
---|
614 | real :: reflon |
---|
615 | real :: scale_top |
---|
616 | real :: xx,yy |
---|
617 | real :: gi2, r2 |
---|
618 | real :: arccos |
---|
619 | |
---|
620 | ! include 'include_wrf_map_utils' |
---|
621 | |
---|
622 | ! begin code |
---|
623 | |
---|
624 | ! compute the reference longitude by rotating 90 degrees to the east |
---|
625 | ! to find the longitude line parallel to the positive x-axis. |
---|
626 | reflon = proj_stdlon + 90. |
---|
627 | |
---|
628 | ! compute numerator term of map scale factor |
---|
629 | scale_top = 1. + proj_hemi * sin(proj_truelat1 * rad_per_deg) |
---|
630 | |
---|
631 | ! compute radius to point of interest |
---|
632 | xx = i - proj_polei |
---|
633 | yy = (j - proj_polej) * proj_hemi |
---|
634 | r2 = xx**2 + yy**2 |
---|
635 | |
---|
636 | ! now the magic code |
---|
637 | if (r2 .eq. 0.) then |
---|
638 | lat = proj_hemi * 90. |
---|
639 | lon = reflon |
---|
640 | else |
---|
641 | gi2 = (proj_rebydx * scale_top)**2. |
---|
642 | lat = deg_per_rad * proj_hemi * asin((gi2-r2)/(gi2+r2)) |
---|
643 | arccos = acos(xx/sqrt(r2)) |
---|
644 | if (yy .gt. 0) then |
---|
645 | lon = reflon + deg_per_rad * arccos |
---|
646 | else |
---|
647 | lon = reflon - deg_per_rad * arccos |
---|
648 | endif |
---|
649 | endif |
---|
650 | |
---|
651 | ! convert to a -180 -> 180 east convention |
---|
652 | if (lon .gt. 180.) lon = lon - 360. |
---|
653 | if (lon .lt. -180.) lon = lon + 360. |
---|
654 | |
---|
655 | return |
---|
656 | end subroutine ijll_ps |
---|
657 | |
---|
658 | |
---|
659 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
660 | subroutine set_lc |
---|
661 | ! initialize the remaining items in the proj structure for a |
---|
662 | ! lambert conformal grid. |
---|
663 | |
---|
664 | use wrf_map_utils_mod |
---|
665 | implicit none |
---|
666 | |
---|
667 | ! include 'include_wrf_map_utils' |
---|
668 | |
---|
669 | real :: arg |
---|
670 | real :: deltalon1 |
---|
671 | real :: tl1r |
---|
672 | real :: ctl1r |
---|
673 | |
---|
674 | ! compute cone factor |
---|
675 | call lc_cone(proj_truelat1, proj_truelat2, proj_cone) |
---|
676 | ! print '(a,f8.6)', 'computed cone factor: ', proj_cone |
---|
677 | ! compute longitude differences and ensure we stay out of the |
---|
678 | ! forbidden "cut zone" |
---|
679 | deltalon1 = proj_lon1 - proj_stdlon |
---|
680 | if (deltalon1 .gt. +180.) deltalon1 = deltalon1 - 360. |
---|
681 | if (deltalon1 .lt. -180.) deltalon1 = deltalon1 + 360. |
---|
682 | |
---|
683 | ! convert truelat1 to radian and compute cos for later use |
---|
684 | tl1r = proj_truelat1 * rad_per_deg |
---|
685 | ctl1r = cos(tl1r) |
---|
686 | |
---|
687 | ! compute the radius to our known lower-left (sw) corner |
---|
688 | proj_rsw = proj_rebydx * ctl1r/proj_cone * & |
---|
689 | (tan((90.*proj_hemi-proj_lat1)*rad_per_deg/2.) / & |
---|
690 | tan((90.*proj_hemi-proj_truelat1)*rad_per_deg/2.)) & |
---|
691 | **proj_cone |
---|
692 | |
---|
693 | ! find pole point |
---|
694 | arg = proj_cone*(deltalon1*rad_per_deg) |
---|
695 | proj_polei = 1. - proj_hemi * proj_rsw * sin(arg) |
---|
696 | proj_polej = 1. + proj_rsw * cos(arg) |
---|
697 | !print '(a,2f10.3)', 'computed pole i/j = ', proj_polei, proj_polej |
---|
698 | |
---|
699 | return |
---|
700 | end subroutine set_lc |
---|
701 | |
---|
702 | |
---|
703 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
704 | subroutine lc_cone(truelat1, truelat2, cone) |
---|
705 | |
---|
706 | ! routine to compute the cone factor of a lambert conformal projection |
---|
707 | |
---|
708 | use wrf_map_utils_mod |
---|
709 | implicit none |
---|
710 | |
---|
711 | ! include 'include_wrf_map_utils' |
---|
712 | |
---|
713 | ! input args |
---|
714 | real :: truelat1 ! (-90 -> 90 degrees n) |
---|
715 | real :: truelat2 ! " " " " " |
---|
716 | |
---|
717 | ! output args |
---|
718 | real :: cone |
---|
719 | |
---|
720 | ! locals |
---|
721 | |
---|
722 | ! begin code |
---|
723 | |
---|
724 | ! first, see if this is a secant or tangent projection. for tangent |
---|
725 | ! projections, truelat1 = truelat2 and the cone is tangent to the |
---|
726 | ! earth surface at this latitude. for secant projections, the cone |
---|
727 | ! intersects the earth surface at each of the distinctly different |
---|
728 | ! latitudes |
---|
729 | if (abs(truelat1-truelat2) .gt. 0.1) then |
---|
730 | |
---|
731 | ! compute cone factor following: |
---|
732 | cone=(alog(cos(truelat1*rad_per_deg)) & |
---|
733 | -alog(cos(truelat2*rad_per_deg))) / & |
---|
734 | (alog(tan((90.-abs(truelat1))*rad_per_deg*0.5 ))- & |
---|
735 | alog(tan((90.-abs(truelat2))*rad_per_deg*0.5 )) ) |
---|
736 | else |
---|
737 | cone = sin(abs(truelat1)*rad_per_deg ) |
---|
738 | endif |
---|
739 | |
---|
740 | return |
---|
741 | end subroutine lc_cone |
---|
742 | |
---|
743 | |
---|
744 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
745 | subroutine ijll_lc( i, j, lat, lon) |
---|
746 | |
---|
747 | ! routine to convert from the (i,j) cartesian coordinate to the |
---|
748 | ! geographical latitude and longitude for a lambert conformal projection. |
---|
749 | |
---|
750 | ! history: |
---|
751 | ! 25 jul 01: corrected by b. shaw, noaa/fsl |
---|
752 | ! |
---|
753 | use wrf_map_utils_mod |
---|
754 | implicit none |
---|
755 | |
---|
756 | ! include 'include_wrf_map_utils' |
---|
757 | |
---|
758 | ! input args |
---|
759 | real :: i ! cartesian x coordinate |
---|
760 | real :: j ! cartesian y coordinate |
---|
761 | |
---|
762 | ! output args |
---|
763 | real :: lat ! latitude (-90->90 deg n) |
---|
764 | real :: lon ! longitude (-180->180 e) |
---|
765 | |
---|
766 | ! locals |
---|
767 | real :: inew |
---|
768 | real :: jnew |
---|
769 | real :: r |
---|
770 | real :: chi,chi1,chi2 |
---|
771 | real :: r2 |
---|
772 | real :: xx |
---|
773 | real :: yy |
---|
774 | |
---|
775 | ! begin code |
---|
776 | |
---|
777 | chi1 = (90. - proj_hemi*proj_truelat1)*rad_per_deg |
---|
778 | chi2 = (90. - proj_hemi*proj_truelat2)*rad_per_deg |
---|
779 | |
---|
780 | ! see if we are in the southern hemispere and flip the indices |
---|
781 | ! if we are. |
---|
782 | if (proj_hemi .eq. -1.) then |
---|
783 | inew = -i + 2. |
---|
784 | jnew = -j + 2. |
---|
785 | else |
---|
786 | inew = i |
---|
787 | jnew = j |
---|
788 | endif |
---|
789 | |
---|
790 | ! compute radius**2 to i/j location |
---|
791 | xx = inew - proj_polei |
---|
792 | yy = proj_polej - jnew |
---|
793 | r2 = (xx*xx + yy*yy) |
---|
794 | r = sqrt(r2)/proj_rebydx |
---|
795 | |
---|
796 | ! convert to lat/lon |
---|
797 | if (r2 .eq. 0.) then |
---|
798 | lat = proj_hemi * 90. |
---|
799 | lon = proj_stdlon |
---|
800 | else |
---|
801 | |
---|
802 | ! longitude |
---|
803 | lon = proj_stdlon + deg_per_rad * atan2(proj_hemi*xx,yy) & |
---|
804 | /proj_cone |
---|
805 | lon = amod(lon+360., 360.) |
---|
806 | |
---|
807 | ! latitude. latitude determined by solving an equation adapted |
---|
808 | ! from: |
---|
809 | ! maling, d.h., 1973: coordinate systems and map projections |
---|
810 | ! equations #20 in appendix i. |
---|
811 | |
---|
812 | if (chi1 .eq. chi2) then |
---|
813 | chi = 2.0*atan( ( r/tan(chi1) )**(1./proj_cone) & |
---|
814 | * tan(chi1*0.5) ) |
---|
815 | else |
---|
816 | chi = 2.0*atan( (r*proj_cone/sin(chi1))**(1./proj_cone) & |
---|
817 | * tan(chi1*0.5)) |
---|
818 | endif |
---|
819 | lat = (90.0-chi*deg_per_rad)*proj_hemi |
---|
820 | |
---|
821 | endif |
---|
822 | |
---|
823 | if (lon .gt. +180.) lon = lon - 360. |
---|
824 | if (lon .lt. -180.) lon = lon + 360. |
---|
825 | |
---|
826 | return |
---|
827 | end subroutine ijll_lc |
---|
828 | |
---|
829 | |
---|
830 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
831 | subroutine llij_lc( lat, lon, i, j) |
---|
832 | |
---|
833 | ! routine to compute the geographical latitude and longitude values |
---|
834 | ! to the cartesian x/y on a lambert conformal projection. |
---|
835 | |
---|
836 | use wrf_map_utils_mod |
---|
837 | implicit none |
---|
838 | ! |
---|
839 | ! include 'include_wrf_map_utils' |
---|
840 | |
---|
841 | ! input args |
---|
842 | real :: lat ! latitude (-90->90 deg n) |
---|
843 | real :: lon ! longitude (-180->180 e) |
---|
844 | |
---|
845 | ! output args |
---|
846 | real :: i ! cartesian x coordinate |
---|
847 | real :: j ! cartesian y coordinate |
---|
848 | |
---|
849 | ! locals |
---|
850 | real :: arg |
---|
851 | real :: deltalon |
---|
852 | real :: tl1r |
---|
853 | real :: rm |
---|
854 | real :: ctl1r |
---|
855 | |
---|
856 | |
---|
857 | ! begin code |
---|
858 | ! print*,proj_clat, proj_clon,proj_stdlon,proj_truelat1,proj_code,proj_cone |
---|
859 | ! compute deltalon between known longitude and standard lon and ensure |
---|
860 | ! it is not in the cut zone |
---|
861 | deltalon = lon - proj_stdlon |
---|
862 | if (deltalon .gt. +180.) deltalon = deltalon - 360. |
---|
863 | if (deltalon .lt. -180.) deltalon = deltalon + 360. |
---|
864 | |
---|
865 | ! convert truelat1 to radian and compute cos for later use |
---|
866 | tl1r = proj_truelat1 * rad_per_deg |
---|
867 | ctl1r = cos(tl1r) |
---|
868 | |
---|
869 | ! radius to desired point |
---|
870 | rm = proj_rebydx * ctl1r/proj_cone * & |
---|
871 | (tan((90.*proj_hemi-lat)*rad_per_deg/2.) / & |
---|
872 | tan((90.*proj_hemi-proj_truelat1)*rad_per_deg/2.))**proj_cone |
---|
873 | |
---|
874 | arg = proj_cone*(deltalon*rad_per_deg) |
---|
875 | i = proj_polei + proj_hemi * rm * sin(arg) |
---|
876 | j = proj_polej - rm * cos(arg) |
---|
877 | |
---|
878 | ! finally, if we are in the southern hemisphere, flip the i/j |
---|
879 | ! values to a coordinate system where (1,1) is the sw corner |
---|
880 | ! (what we assume) which is different than the original ncep |
---|
881 | ! algorithms which used the ne corner as the origin in the |
---|
882 | ! southern hemisphere (left-hand vs. right-hand coordinate?) |
---|
883 | if (proj_hemi .eq. -1.) then |
---|
884 | i = 2. - i |
---|
885 | j = 2. - j |
---|
886 | endif |
---|
887 | |
---|
888 | return |
---|
889 | end subroutine llij_lc |
---|
890 | |
---|
891 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
892 | SUBROUTINE set_merc |
---|
893 | |
---|
894 | ! Sets up the remaining basic elements for the mercator projection |
---|
895 | use wrf_map_utils_mod |
---|
896 | |
---|
897 | IMPLICIT NONE |
---|
898 | ! TYPE(proj_info), INTENT(INOUT) :: proj |
---|
899 | REAL :: clain |
---|
900 | |
---|
901 | ! Preliminary variables |
---|
902 | |
---|
903 | clain = COS(rad_per_deg*proj_truelat1) |
---|
904 | proj_dlon = proj_dx / (earth_radius_m * clain) |
---|
905 | |
---|
906 | ! Compute distance from equator to origin, and store in the |
---|
907 | ! proj%rsw tag. |
---|
908 | |
---|
909 | proj_rsw = 0. |
---|
910 | IF (proj_lat1 .NE. 0.) THEN |
---|
911 | proj_rsw = (ALOG(TAN(0.5*((proj_lat1+90.)*rad_per_deg))))/proj_dlon |
---|
912 | ENDIF |
---|
913 | RETURN |
---|
914 | END SUBROUTINE set_merc |
---|
915 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
916 | |
---|
917 | SUBROUTINE llij_merc(lat, lon, i, j) |
---|
918 | |
---|
919 | ! Compute i/j coordinate from lat lon for mercator projection |
---|
920 | use wrf_map_utils_mod |
---|
921 | |
---|
922 | IMPLICIT NONE |
---|
923 | REAL, INTENT(IN) :: lat |
---|
924 | REAL, INTENT(IN) :: lon |
---|
925 | ! TYPE(proj_info),INTENT(IN) :: proj |
---|
926 | REAL,INTENT(OUT) :: i |
---|
927 | REAL,INTENT(OUT) :: j |
---|
928 | REAL :: deltalon |
---|
929 | |
---|
930 | deltalon = lon - proj_lon1 |
---|
931 | ! deltalon = lon - proj_stdlon |
---|
932 | ! JB modif |
---|
933 | ! IF (deltalon .LT. -180.) deltalon = deltalon + 360. |
---|
934 | ! IF (deltalon .GT. 180.) deltalon = deltalon - 360. |
---|
935 | i = 1. + (deltalon/(proj_dlon*deg_per_rad)) |
---|
936 | j = 1. + (ALOG(TAN(0.5*((lat + 90.) * rad_per_deg)))) / & |
---|
937 | proj_dlon - proj_rsw |
---|
938 | ! print*,'IN MERCATOR',lat,lon,i,j,proj_lon1,deltalon |
---|
939 | ! pause |
---|
940 | RETURN |
---|
941 | END SUBROUTINE llij_merc |
---|
942 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
943 | SUBROUTINE ijll_merc(i, j, lat, lon) |
---|
944 | |
---|
945 | ! Compute the lat/lon from i/j for mercator projection |
---|
946 | use wrf_map_utils_mod |
---|
947 | |
---|
948 | IMPLICIT NONE |
---|
949 | REAL,INTENT(IN) :: i |
---|
950 | REAL,INTENT(IN) :: j |
---|
951 | ! TYPE(proj_info),INTENT(IN) :: proj |
---|
952 | REAL, INTENT(OUT) :: lat |
---|
953 | REAL, INTENT(OUT) :: lon |
---|
954 | |
---|
955 | |
---|
956 | lat = 2.0*ATAN(EXP(proj_dlon*(proj_rsw + j-1.)))*deg_per_rad - 90. |
---|
957 | lon = (i-1.)*proj_dlon*deg_per_rad + proj_lon1 |
---|
958 | IF (lon.GT.180.) lon = lon - 360. |
---|
959 | IF (lon.LT.-180.) lon = lon + 360. |
---|
960 | RETURN |
---|
961 | END SUBROUTINE ijll_merc |
---|
962 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
963 | SUBROUTINE llij_latlon(lat, lon, i, j) |
---|
964 | |
---|
965 | ! Compute the i/j location of a lat/lon on a LATLON grid. |
---|
966 | use wrf_map_utils_mod |
---|
967 | |
---|
968 | IMPLICIT NONE |
---|
969 | REAL, INTENT(IN) :: lat |
---|
970 | REAL, INTENT(IN) :: lon |
---|
971 | ! TYPE(proj_info), INTENT(IN) :: proj |
---|
972 | REAL, INTENT(OUT) :: i |
---|
973 | REAL, INTENT(OUT) :: j |
---|
974 | |
---|
975 | REAL :: deltalat |
---|
976 | REAL :: deltalon |
---|
977 | REAL :: lon360 |
---|
978 | REAL :: latinc |
---|
979 | REAL :: loninc |
---|
980 | |
---|
981 | ! Extract the latitude and longitude increments for this grid |
---|
982 | ! (e.g., 2.5 deg for NCEP reanalysis) from the proj structure, where |
---|
983 | ! loninc is saved in the stdlon tag and latinc is saved in truelat1 |
---|
984 | |
---|
985 | latinc = proj_truelat1 |
---|
986 | loninc = proj_stdlon |
---|
987 | |
---|
988 | ! Compute deltalat and deltalon as the difference between the input |
---|
989 | ! lat/lon and the origin lat/lon |
---|
990 | |
---|
991 | deltalat = lat - proj_lat1 |
---|
992 | ! To account for issues around the dateline, convert the incoming |
---|
993 | ! longitudes to be 0->360. |
---|
994 | IF (lon .LT. 0) THEN |
---|
995 | lon360 = lon + 360. |
---|
996 | ELSE |
---|
997 | lon360 = lon |
---|
998 | ENDIF |
---|
999 | deltalon = lon360 - proj_lon1 |
---|
1000 | IF (deltalon .LT. 0) deltalon = deltalon + 360. |
---|
1001 | |
---|
1002 | ! Compute i/j |
---|
1003 | i = deltalon/loninc + 1. |
---|
1004 | j = deltalat/latinc + 1. |
---|
1005 | RETURN |
---|
1006 | END SUBROUTINE llij_latlon |
---|
1007 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
1008 | SUBROUTINE ijll_latlon(i, j, lat, lon) |
---|
1009 | |
---|
1010 | ! Compute the lat/lon location of an i/j on a LATLON grid. |
---|
1011 | use wrf_map_utils_mod |
---|
1012 | |
---|
1013 | IMPLICIT NONE |
---|
1014 | REAL, INTENT(IN) :: i |
---|
1015 | REAL, INTENT(IN) :: j |
---|
1016 | ! TYPE(proj_info), INTENT(IN) :: proj |
---|
1017 | REAL, INTENT(OUT) :: lat |
---|
1018 | REAL, INTENT(OUT) :: lon |
---|
1019 | |
---|
1020 | REAL :: deltalat |
---|
1021 | REAL :: deltalon |
---|
1022 | REAL :: lon360 |
---|
1023 | REAL :: latinc |
---|
1024 | REAL :: loninc |
---|
1025 | |
---|
1026 | ! Extract the latitude and longitude increments for this grid |
---|
1027 | ! (e.g., 2.5 deg for NCEP reanalysis) from the proj structure, where |
---|
1028 | ! loninc is saved in the stdlon tag and latinc is saved in truelat1 |
---|
1029 | |
---|
1030 | latinc = proj_truelat1 |
---|
1031 | loninc = proj_stdlon |
---|
1032 | |
---|
1033 | ! Compute deltalat and deltalon |
---|
1034 | |
---|
1035 | deltalat = (j-1.)*latinc |
---|
1036 | deltalon = (i-1.)*loninc |
---|
1037 | lat = proj_lat1 + deltalat |
---|
1038 | lon = proj_lon1 + deltalon |
---|
1039 | |
---|
1040 | IF ((ABS(lat) .GT. 90.).OR.(ABS(deltalon) .GT.360.)) THEN |
---|
1041 | ! Off the earth for this grid |
---|
1042 | lat = -999. |
---|
1043 | lon = -999. |
---|
1044 | ELSE |
---|
1045 | lon = lon + 360. |
---|
1046 | lon = AMOD(lon,360.) |
---|
1047 | IF (lon .GT. 180.) lon = lon -360. |
---|
1048 | ENDIF |
---|
1049 | |
---|
1050 | RETURN |
---|
1051 | END SUBROUTINE ijll_latlon |
---|
1052 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
1053 | SUBROUTINE gridwind_to_truewind(lon,ugrid,vgrid,utrue,vtrue) |
---|
1054 | |
---|
1055 | ! Subroutine to convert a wind from grid north to true north. |
---|
1056 | use wrf_map_utils_mod |
---|
1057 | |
---|
1058 | IMPLICIT NONE |
---|
1059 | |
---|
1060 | ! Arguments |
---|
1061 | REAL, INTENT(IN) :: lon ! Longitude of point in degrees |
---|
1062 | ! TYPE(proj_info),INTENT(IN):: proj ! Projection info structure |
---|
1063 | REAL, INTENT(IN) :: ugrid ! U-component, grid-relative |
---|
1064 | REAL, INTENT(IN) :: vgrid ! V-component, grid-relative |
---|
1065 | REAL, INTENT(OUT) :: utrue ! U-component, earth-relative |
---|
1066 | REAL, INTENT(OUT) :: vtrue ! V-component, earth-relative |
---|
1067 | |
---|
1068 | ! Locals |
---|
1069 | REAL :: alpha |
---|
1070 | REAL :: diff |
---|
1071 | |
---|
1072 | IF ((proj_code .EQ. PROJ_PS).OR.(proj_code .EQ. PROJ_LC))THEN |
---|
1073 | diff = lon - proj_stdlon |
---|
1074 | IF (diff .GT. 180.) diff = diff - 360. |
---|
1075 | IF (diff .LT.-180.) diff = diff + 360. |
---|
1076 | alpha = diff * proj_cone * rad_per_deg * SIGN(1.,proj_truelat1) |
---|
1077 | utrue = vgrid * SIN(alpha) + ugrid * COS(alpha) |
---|
1078 | vtrue = vgrid * COS(alpha) - ugrid * SIN(alpha) |
---|
1079 | ELSEIF ((proj_code .EQ. PROJ_MERC).OR.(proj_code .EQ.PROJ_LATLON))THEN |
---|
1080 | utrue = ugrid |
---|
1081 | vtrue = vgrid |
---|
1082 | ELSE |
---|
1083 | PRINT '(A)', 'Unrecognized projection.' |
---|
1084 | STOP 'GRIDWIND_TO_TRUEWIND' |
---|
1085 | ENDIF |
---|
1086 | |
---|
1087 | RETURN |
---|
1088 | END SUBROUTINE gridwind_to_truewind |
---|
1089 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
1090 | |
---|