[16] | 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 | |
---|