| 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 | |
|---|
