geogrid/src/llxy_module.F

   1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   2 ! MODULE LLXY_MODULE
   3 !
   4 ! This module handles transformations between model grid coordinates and
   5 !   latitude-longitude coordinates. The actual transformations are done through
   6 !   the map_utils module.
   7 !
   8 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   9 module llxy_module
  10
  11    use gridinfo_module
  12    use list_module
  13    use map_utils
  14    use module_debug
  15    use misc_definitions_module
  16
  17    ! Parameters
  18    integer, parameter :: MAX_SOURCE_LEVELS = 20
  19
  20    ! Variables
  21    integer :: current_nest_number
  22    integer :: SOURCE_PROJ = 0
  23    ! The following arrays hold values for all available domains
  24    ! NOTE: The entries in the arrays for "domain 0" are used for projection
  25    !       information of user-specified source data
  26    type (proj_info), dimension(-MAX_SOURCE_LEVELS:MAX_DOMAINS) :: proj_stack
  27
  28    ! The projection and domain that we have computed constants for
  29    integer :: computed_proj = INVALID
  30    integer :: computed_domain = INVALID
  31
  32    contains
  33
  34    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  35    ! Name: push_source_projection
  36    !
  37    ! Purpose:
  38    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  39    subroutine push_source_projection(iprojection, user_stand_lon, user_truelat1, user_truelat2, &
  40                                   user_dxkm, user_dykm, user_dlat, user_dlon, user_known_x, &
  41                                   user_known_y, user_known_lat, user_known_lon, &
  42                                   user_pole_lat, user_pole_lon, &
  43                                   user_centerlat, user_centerlon, &
  44                                   user_centeri, user_centerj, &
  45                                   earth_radius)
  46
  47       implicit none
  48
  49       ! Arguments
  50       integer, intent(in) :: iprojection
  51       real, intent(in) :: user_stand_lon, user_truelat1, user_truelat2, user_dxkm, user_dykm, &
  52                           user_dlat, user_dlon, &
  53                           user_known_x, user_known_y, user_known_lat, user_known_lon
  54       real, intent(in), optional :: earth_radius
  55       real, intent(in), optional :: user_centerlon, user_centerlat, user_pole_lat, user_pole_lon
  56       real, intent(in), optional :: user_centerj, user_centeri
  57
  58
  59       SOURCE_PROJ = SOURCE_PROJ-1
  60       if (SOURCE_PROJ < -MAX_SOURCE_LEVELS) then
  61          call mprintf(.true.,ERROR,'In push_user_projection(), too many levels of user projections.')
  62       end if
  63
  64       call map_init(proj_stack(SOURCE_PROJ))
  65
  66       if (iprojection == PROJ_LATLON) then
  67          call map_set(iprojection, proj_stack(SOURCE_PROJ), &
  68                       lat1=user_known_lat, &
  69                       lon1=user_known_lon, &
  70                       knowni=user_known_x, &
  71                       knownj=user_known_y, &
  72                       nxmax=nint(360.0 / user_dlon), &
  73                       latinc=user_dlat, &
  74                       loninc=user_dlon, &
  75                       r_earth=earth_radius)
  76
  77       else if (iprojection == PROJ_MERC) then
  78          call map_set(iprojection, proj_stack(SOURCE_PROJ), &
  79                       truelat1=user_truelat1, &
  80                       lat1=user_known_lat, &
  81                       lon1=user_known_lon, &
  82                       knowni=user_known_x, &
  83                       knownj=user_known_y, &
  84                       dx=user_dxkm, &
  85                       r_earth=earth_radius)
  86
  87       else if (iprojection == PROJ_CYL) then
  88          call mprintf(.true.,ERROR,'Should not have PROJ_CYL as projection for ' &
  89                           //'source data in push_source_projection()')
  90
  91       else if (iprojection == PROJ_CASSINI) then
  92
  93          call map_set(iprojection, proj_stack(SOURCE_PROJ), &
  94                       latinc=user_dlat, &
  95                       loninc=user_dlon, &
  96                       stdlon=user_stand_lon, &
  97                       lat1=user_centerlat, &
  98                       lon1=user_centerlon, &
  99                       lat0=user_pole_lat, &
 100                       lon0=user_pole_lon, &
 101                       knowni=user_centeri, &
 102                       knownj=user_centerj, &
 103                       r_earth=earth_radius)
 104
 105       else if (iprojection == PROJ_LC) then
 106          call map_set(iprojection, proj_stack(SOURCE_PROJ), &
 107                       truelat1=user_truelat1, &
 108                       truelat2=user_truelat2, &
 109                       stdlon=user_stand_lon, &
 110                       lat1=user_known_lat, &
 111                       lon1=user_known_lon, &
 112                       knowni=user_known_x, &
 113                       knownj=user_known_y, &
 114                       dx=user_dxkm, &
 115                       r_earth=earth_radius)
 116
 117       else if (iprojection == PROJ_ALBERS_NAD83) then
 118          call map_set(iprojection, proj_stack(SOURCE_PROJ), &
 119                       truelat1=user_truelat1, &
 120                       truelat2=user_truelat2, &
 121                       stdlon=user_stand_lon, &
 122                       lat1=user_known_lat, &
 123                       lon1=user_known_lon, &
 124                       knowni=user_known_x, &
 125                       knownj=user_known_y, &
 126                       dx=user_dxkm, &
 127                       r_earth=earth_radius)
 128
 129       else if (iprojection == PROJ_PS) then
 130          call map_set(iprojection, proj_stack(SOURCE_PROJ), &
 131                       truelat1=user_truelat1, &
 132                       stdlon=user_stand_lon, &
 133                       lat1=user_known_lat, &
 134                       lon1=user_known_lon, &
 135                       knowni=user_known_x, &
 136                       knownj=user_known_y, &
 137                       dx=user_dxkm, &
 138                       r_earth=earth_radius)
 139
 140       else if (iprojection == PROJ_PS_WGS84) then
 141          call map_set(iprojection, proj_stack(SOURCE_PROJ), &
 142                       truelat1=user_truelat1, &
 143                       stdlon=user_stand_lon, &
 144                       lat1=user_known_lat, &
 145                       lon1=user_known_lon, &
 146                       knowni=user_known_x, &
 147                       knownj=user_known_y, &
 148                       dx=user_dxkm, &
 149                       r_earth=earth_radius)
 150
 151       else if (iprojection == PROJ_GAUSS) then
 152          call map_set(iprojection, proj_stack(SOURCE_PROJ), &
 153                       lat1=user_known_lat, &
 154                       lon1=user_known_lon, &
 155                       nxmax=nint(360.0 / user_dlon), &
 156                       nlat=nint(user_dlat), &
 157                       loninc=user_dlon, &
 158                       r_earth=earth_radius)
 159
 160       else if (iprojection == PROJ_ROTLL) then
 161   ! BUG: Implement this projection.
 162
 163       end if
 164
 165    end subroutine push_source_projection
 166
 167
 168    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 169    ! Name: pop_source_projection
 170    !
 171    ! Purpose:
 172    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 173    subroutine pop_source_projection()
 174
 175       implicit none
 176
 177       SOURCE_PROJ = SOURCE_PROJ+1
 178
 179       call mprintf((SOURCE_PROJ > 0), ERROR, &
 180                    'In pop_user_projection(), projection stack has overflowed.')
 181
 182    end subroutine pop_source_projection
 183
 184
 185 #ifdef _METGRID
 186    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 187    ! Name: set_domain_projection
 188    !
 189    ! Purpose:
 190    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 191    subroutine set_domain_projection(iprojection, user_stand_lon, user_truelat1, user_truelat2, &
 192                                   user_dxkm, user_dykm, user_dlat, user_dlon, &
 193                                   user_xdim, user_ydim, user_known_x, &
 194                                   user_known_y, user_known_lat, user_known_lon, &
 195                                   user_pole_lat, user_pole_lon, earth_radius)
 196
 197       implicit none
 198
 199       ! Arguments
 200       integer, intent(in) :: iprojection
 201       integer, intent(in) :: user_xdim, user_ydim
 202       real, intent(in) :: user_stand_lon, user_truelat1, user_truelat2, &
 203                           user_dxkm, user_dykm, user_dlat, user_dlon, &
 204                           user_known_x, user_known_y, user_known_lat, user_known_lon, &
 205                           user_pole_lat, user_pole_lon
 206       real, intent(in), optional :: earth_radius
 207
 208       current_nest_number = 1
 209
 210       call map_init(proj_stack(current_nest_number))
 211
 212       if (iprojection == PROJ_LATLON) then
 213          call map_set(iprojection, proj_stack(current_nest_number), &
 214                       lat1=user_known_lat, &
 215                       lon1=user_known_lon, &
 216                       latinc=user_dlat, &
 217                       loninc=user_dlon, &
 218                       r_earth=earth_radius)
 219
 220       else if (iprojection == PROJ_MERC) then
 221          call map_set(iprojection, proj_stack(current_nest_number), &
 222                       truelat1=user_truelat1, &
 223                       lat1=user_known_lat, &
 224                       lon1=user_known_lon, &
 225                       knowni=user_known_x, &
 226                       knownj=user_known_y, &
 227                       dx=user_dxkm, &
 228                       r_earth=earth_radius)
 229
 230       else if (iprojection == PROJ_CYL) then
 231          call map_set(iprojection, proj_stack(current_nest_number), &
 232                       latinc=user_dlat, &
 233                       loninc=user_dlon, &
 234                       stdlon=user_stand_lon, &
 235                       r_earth=earth_radius)
 236
 237       else if (iprojection == PROJ_CASSINI) then
 238          call map_set(iprojection, proj_stack(current_nest_number), &
 239                       latinc=user_dlat, &
 240                       loninc=user_dlon, &
 241                       dx=user_dxkm,        &
 242                       dy=user_dykm,        &
 243                       stdlon=user_stand_lon, &
 244                       lat1=user_known_lat, &
 245                       lon1=user_known_lon, &
 246                       lat0=user_pole_lat, &
 247                       lon0=user_pole_lon, &
 248                       knowni=user_known_x, &
 249                       knownj=user_known_y, &
 250                       r_earth=earth_radius)
 251
 252       else if (iprojection == PROJ_LC) then
 253          call map_set(iprojection, proj_stack(current_nest_number), &
 254                       truelat1=user_truelat1, &
 255                       truelat2=user_truelat2, &
 256                       stdlon=user_stand_lon, &
 257                       lat1=user_known_lat, &
 258                       lon1=user_known_lon, &
 259                       knowni=user_known_x, &
 260                       knownj=user_known_y, &
 261                       dx=user_dxkm, &
 262                       r_earth=earth_radius)
 263
 264       else if (iprojection == PROJ_ALBERS_NAD83) then
 265          call map_set(iprojection, proj_stack(current_nest_number), &
 266                       truelat1=user_truelat1, &
 267                       truelat2=user_truelat2, &
 268                       stdlon=user_stand_lon, &
 269                       lat1=user_known_lat, &
 270                       lon1=user_known_lon, &
 271                       knowni=user_known_x, &
 272                       knownj=user_known_y, &
 273                       dx=user_dxkm, &
 274                       r_earth=earth_radius)
 275
 276       else if (iprojection == PROJ_PS) then
 277          call map_set(iprojection, proj_stack(current_nest_number), &
 278                       truelat1=user_truelat1, &
 279                       stdlon=user_stand_lon, &
 280                       lat1=user_known_lat, &
 281                       lon1=user_known_lon, &
 282                       knowni=user_known_x, &
 283                       knownj=user_known_y, &
 284                       dx=user_dxkm, &
 285                       r_earth=earth_radius)
 286
 287       else if (iprojection == PROJ_PS_WGS84) then
 288          call map_set(iprojection, proj_stack(current_nest_number), &
 289                       truelat1=user_truelat1, &
 290                       stdlon=user_stand_lon, &
 291                       lat1=user_known_lat, &
 292                       lon1=user_known_lon, &
 293                       knowni=user_known_x, &
 294                       knownj=user_known_y, &
 295                       dx=user_dxkm)
 296
 297       else if (iprojection == PROJ_GAUSS) then
 298          call map_set(iprojection, proj_stack(current_nest_number), &
 299                       lat1=user_known_lat, &
 300                       lon1=user_known_lon, &
 301                       nlat=nint(user_dlat), &
 302                       loninc=user_dlon, &
 303                       r_earth=earth_radius)
 304
 305       else if (iprojection == PROJ_ROTLL) then
 306          call map_set(iprojection, proj_stack(current_nest_number), &
 307                       ixdim=user_xdim, &
 308                       jydim=user_ydim, &
 309                       phi=user_dlat, &
 310                       lambda=user_dlon, &
 311                       lat1=user_known_lat, &
 312                       lon1=user_known_lon, &
 313                       stagger=HH, &
 314                       latinc=user_dykm, &
 315                       loninc=user_dxkm, &
 316                       r_earth=earth_radius)
 317
 318       end if
 319
 320    end subroutine set_domain_projection
 321 #endif
 322
 323
 324 #ifdef _GEOGRID
 325    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 326    ! Name: compute_nest_locations
 327    !
 328    ! Purpose: This routine computes the variables necessary in determining the
 329    !   location of all nests without reference to the parent or coarse domains.
 330    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 331    subroutine compute_nest_locations()
 332
 333       implicit none
 334
 335       ! Local variables
 336       integer :: i
 337       real :: temp_known_x, temp_known_y, temp_known_lat, temp_known_lon, &
 338               temp_dxkm, temp_dykm, temp_dlat, temp_dlon
 339
 340       ! Set location of coarse/mother domain
 341       call map_init(proj_stack(1))
 342
 343       if (iproj_type == PROJ_LATLON) then
 344          call map_set(iproj_type, proj_stack(1), &
 345                       lat1=known_lat, &
 346                       lon1=known_lon, &
 347                       latinc=dykm, &
 348                       loninc=dxkm)
 349
 350       else if (iproj_type == PROJ_MERC) then
 351          call map_set(iproj_type, proj_stack(1), &
 352                       truelat1=truelat1, &
 353                       lat1=known_lat, &
 354                       lon1=known_lon, &
 355                       knowni=known_x, &
 356                       knownj=known_y, &
 357                       dx=dxkm)
 358
 359       else if (iproj_type == PROJ_CYL) then
 360          call map_set(iproj_type, proj_stack(1), &
 361                       latinc=dlatdeg, &
 362                       loninc=dlondeg, &
 363                       stdlon=stand_lon)
 364
 365       else if (iproj_type == PROJ_CASSINI) then
 366          call map_set(iproj_type, proj_stack(1), &
 367                       latinc=dlatdeg, &
 368                       loninc=dlondeg, &
 369                       dx=dxkm,       &
 370                       dy=dykm,       &
 371                       stdlon=stand_lon, &
 372                       knowni=known_x, &
 373                       knownj=known_y, &
 374                       lat0=pole_lat, &
 375                       lon0=pole_lon, &
 376                       lat1=known_lat, &
 377                       lon1=known_lon)
 378
 379       else if (iproj_type == PROJ_LC) then
 380          call map_set(iproj_type, proj_stack(1), &
 381                       truelat1=truelat1, &
 382                       truelat2=truelat2, &
 383                       stdlon=stand_lon, &
 384                       lat1=known_lat, &
 385                       lon1=known_lon, &
 386                       knowni=known_x, &
 387                       knownj=known_y, &
 388                       dx=dxkm)
 389
 390       else if (iproj_type == PROJ_ALBERS_NAD83) then
 391          call map_set(iproj_type, proj_stack(1), &
 392                       truelat1=truelat1, &
 393                       truelat2=truelat2, &
 394                       stdlon=stand_lon, &
 395                       lat1=known_lat, &
 396                       lon1=known_lon, &
 397                       knowni=known_x, &
 398                       knownj=known_y, &
 399                       dx=dxkm)
 400
 401       else if (iproj_type == PROJ_PS) then
 402          call map_set(iproj_type, proj_stack(1), &
 403                       truelat1=truelat1, &
 404                       stdlon=stand_lon, &
 405                       lat1=known_lat, &
 406                       lon1=known_lon, &
 407                       knowni=known_x, &
 408                       knownj=known_y, &
 409                       dx=dxkm)
 410
 411       else if (iproj_type == PROJ_PS_WGS84) then
 412          call map_set(iproj_type, proj_stack(1), &
 413                       truelat1=truelat1, &
 414                       stdlon=stand_lon, &
 415                       lat1=known_lat, &
 416                       lon1=known_lon, &
 417                       knowni=known_x, &
 418                       knownj=known_y, &
 419                       dx=dxkm)
 420
 421       else if (iproj_type == PROJ_GAUSS) then
 422          call map_set(iproj_type, proj_stack(current_nest_number), &
 423                       lat1=known_lat, &
 424                       lon1=known_lon, &
 425                       nlat=nint(dykm), &
 426                       loninc=dxkm)
 427
 428       else if (iproj_type == PROJ_ROTLL) then
 429          call map_set(iproj_type, proj_stack(1), &
 430                       ixdim=ixdim(1), &
 431                       jydim=jydim(1), &
 432                       phi=phi, &
 433                       lambda=lambda, &
 434                       lat1=known_lat, &
 435                       lon1=known_lon, &
 436                       latinc=dykm, &
 437                       loninc=dxkm, &
 438                       stagger=HH)
 439
 440       end if
 441
 442       ! Now we can compute lat/lon <-> x/y for coarse domain
 443       call select_domain(1)
 444
 445       ! Call a recursive procedure to find the lat/lon of the centerpoint for
 446       !   each domain
 447       do i=2,n_domains
 448
 449          temp_known_x = real(ixdim(i))/2.
 450          temp_known_y = real(jydim(i))/2.
 451
 452          call find_known_latlon(i, temp_known_x, temp_known_y, &
 453                                 temp_known_lat, temp_known_lon, &
 454                                 temp_dxkm, temp_dykm, temp_dlat, temp_dlon)
 455
 456          if (iproj_type == PROJ_LATLON) then
 457             call map_set(iproj_type, proj_stack(i), &
 458                          lat1=temp_known_lat, &
 459                          lon1=temp_known_lon, &
 460                          latinc=temp_dlat, &
 461                          loninc=temp_dlon)
 462
 463          else if (iproj_type == PROJ_MERC) then
 464             call map_set(iproj_type, proj_stack(i), &
 465                          truelat1=truelat1, &
 466                          lat1=temp_known_lat, &
 467                          lon1=temp_known_lon, &
 468                          knowni=temp_known_x, &
 469                          knownj=temp_known_y, &
 470                          dx=temp_dxkm)
 471
 472          else if (iproj_type == PROJ_CYL) then
 473             call mprintf(.true.,ERROR,'Don''t know how to do nesting with PROJ_CYL ' &
 474                                       //'in compute_nest_locations()')
 475
 476          else if (iproj_type == PROJ_CASSINI) then
 477             call map_set(iproj_type, proj_stack(i), &
 478                          latinc=temp_dlat, &
 479                          loninc=temp_dlon, &
 480                          dx=temp_dxkm,  &
 481                          dy=temp_dykm,  &
 482                          stdlon=stand_lon, &
 483                          knowni=temp_known_x, &
 484                          knownj=temp_known_y, &
 485                          lat0=pole_lat, &
 486                          lon0=pole_lon, &
 487                          lat1=temp_known_lat, &
 488                          lon1=temp_known_lon)
 489
 490          else if (iproj_type == PROJ_LC) then
 491             call map_set(iproj_type, proj_stack(i), &
 492                          truelat1=truelat1, &
 493                          truelat2=truelat2, &
 494                          stdlon=stand_lon, &
 495                          lat1=temp_known_lat, &
 496                          lon1=temp_known_lon, &
 497                          knowni=temp_known_x, &
 498                          knownj=temp_known_y, &
 499                          dx=temp_dxkm)
 500
 501          else if (iproj_type == PROJ_ALBERS_NAD83) then
 502             call map_set(iproj_type, proj_stack(i), &
 503                          truelat1=truelat1, &
 504                          truelat2=truelat2, &
 505                          stdlon=stand_lon, &
 506                          lat1=temp_known_lat, &
 507                          lon1=temp_known_lon, &
 508                          knowni=temp_known_x, &
 509                          knownj=temp_known_y, &
 510                          dx=temp_dxkm)
 511
 512          else if (iproj_type == PROJ_PS) then
 513             call map_set(iproj_type, proj_stack(i), &
 514                          truelat1=truelat1, &
 515                          stdlon=stand_lon, &
 516                          lat1=temp_known_lat, &
 517                          lon1=temp_known_lon, &
 518                          knowni=temp_known_x, &
 519                          knownj=temp_known_y, &
 520                          dx=temp_dxkm)
 521
 522          else if (iproj_type == PROJ_PS_WGS84) then
 523             call map_set(iproj_type, proj_stack(i), &
 524                          truelat1=truelat1, &
 525                          stdlon=stand_lon, &
 526                          lat1=temp_known_lat, &
 527                          lon1=temp_known_lon, &
 528                          knowni=temp_known_x, &
 529                          knownj=temp_known_y, &
 530                          dx=temp_dxkm)
 531
 532          else if (iproj_type == PROJ_GAUSS) then
 533             call map_set(iproj_type, proj_stack(current_nest_number), &
 534                          lat1=temp_known_lat, &
 535                          lon1=temp_known_lon, &
 536                          nlat=nint(temp_dykm), &
 537                          loninc=temp_dxkm)
 538
 539          else if (iproj_type == PROJ_ROTLL) then
 540    ! BUG: Implement this projection.
 541
 542          end if
 543
 544       end do
 545
 546    end subroutine compute_nest_locations
 547
 548
 549    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 550    ! Name: find_known_latlon
 551    !
 552    ! Purpose: This recursive routine computes the latitude and longitude for a
 553    !   specified x/y location in the given nest number, and also computes the
 554    !   grid spacing
 555    !
 556    ! NOTE: This routine assumes that xytoll will work correctly for the
 557    !       coarse domain.
 558    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 559    recursive subroutine find_known_latlon(n, rx, ry, rlat, rlon, dx, dy, dlat, dlon)
 560
 561       implicit none
 562
 563       ! Arguments
 564       integer, intent(in) :: n
 565       real, intent(in) :: rx, ry
 566       real, intent(out) :: rlat, rlon, dx, dy, dlat, dlon
 567
 568       ! Local variables
 569       real :: x_in_parent, y_in_parent
 570
 571       if (n == 1) then   ! Stopping case for the recursion
 572
 573          dx = dxkm
 574          dy = dykm
 575          dlat = dlatdeg
 576          dlon = dlondeg
 577          call ij_to_latlon(proj_stack(current_nest_number), rx, ry, rlat, rlon)
 578
 579          return
 580
 581       else               ! Recursive case
 582
 583          x_in_parent = (rx - ((parent_grid_ratio(n)+1.)/2.)) &
 584                       / parent_grid_ratio(n) + parent_ll_x(n)
 585          y_in_parent = (ry - ((parent_grid_ratio(n)+1.)/2.)) &
 586                       / parent_grid_ratio(n) + parent_ll_y(n)
 587
 588          call find_known_latlon(parent_id(n), x_in_parent, y_in_parent, rlat, rlon, dx, dy, dlat, dlon)
 589
 590          dx = dx / parent_grid_ratio(n)
 591          dy = dy / parent_grid_ratio(n)
 592          dlat = dlat / parent_grid_ratio(n)
 593          dlon = dlon / parent_grid_ratio(n)
 594       end if
 595
 596    end subroutine find_known_latlon
 597
 598
 599    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 600    ! Name: compute_nest_level_info
 601    !
 602    ! Purpose: This routine computes the parameters describing a nesting level for
 603    !          NMM grids.
 604    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 605    subroutine compute_nest_level_info()
 606
 607       implicit none
 608
 609       ! Local variables
 610       integer :: i, nest_level, temp
 611       type (list) :: level_list
 612
 613       call list_init(level_list)
 614
 615       ! Set location of coarse/mother domain
 616       call map_init(proj_stack(1))
 617
 618       call map_set(PROJ_ROTLL, proj_stack(1), &
 619                    ixdim=ixdim(1), &
 620                    jydim=jydim(1), &
 621                    phi=phi, &
 622                    lambda=lambda, &
 623                    lat1=known_lat, &
 624                    lon1=known_lon, &
 625                    latinc=dykm, &
 626                    loninc=dxkm, &
 627                    stagger=HH)
 628
 629       parent_ur_x(1) = real(ixdim(1))
 630       parent_ur_y(1) = real(jydim(1))
 631
 632       do i=2,n_domains
 633
 634          nest_level = get_nest_level(i)
 635
 636          if (.not. list_search(level_list, ikey=nest_level, ivalue=temp)) then
 637
 638             call list_insert(level_list, ikey=nest_level, ivalue=nest_level)
 639
 640             ixdim(nest_level) = ixdim(1)*(3**(nest_level-1))-(3**(nest_level-1)-1)
 641             jydim(nest_level) = jydim(1)*(3**(nest_level-1))-(3**(nest_level-1)-1)
 642
 643             parent_ur_x(nest_level) = ixdim(nest_level)
 644             parent_ur_y(nest_level) = jydim(nest_level)
 645
 646             call map_set(PROJ_ROTLL, proj_stack(nest_level), &
 647                          ixdim = ixdim(nest_level), &
 648                          jydim = jydim(nest_level), &
 649                          phi    = phi, &
 650                          lambda = lambda, &
 651                          lat1=known_lat, &
 652                          lon1=known_lon, &
 653                          latinc=(dykm/real((3**(nest_level-1)))), &
 654                          loninc=(dxkm/real((3**(nest_level-1)))), &
 655                          stagger=HH)
 656          end if
 657
 658       end do
 659
 660       call list_destroy(level_list)
 661
 662    end subroutine compute_nest_level_info
 663
 664
 665    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 666    ! Name: get_domain_resolution
 667    !
 668    ! Purpose:
 669    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 670    subroutine get_domain_resolution(dom_dx, dom_dy)
 671
 672       implicit none
 673
 674       ! Arguments
 675       real, intent(out) :: dom_dx, dom_dy
 676
 677       ! The proj_info structure only stores dx, so set both dom_dx and dom_dy to dx
 678       dom_dx = proj_stack(current_nest_number)%dx
 679       dom_dy = proj_stack(current_nest_number)%dx
 680
 681    end subroutine get_domain_resolution
 682
 683
 684    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 685    ! Name: get_nest_level
 686    !
 687    ! Purpose: This function returns, given a grid ID number, the nesting level of
 688    !   that domain; the coarse domain is taken to have nesting level 1.
 689    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 690    function get_nest_level(i)
 691
 692       implicit none
 693
 694       ! Arguments
 695       integer, intent(in) :: i
 696
 697       ! Local variables
 698       integer :: j
 699
 700       ! Return value
 701       integer :: get_nest_level
 702
 703       ! If argument is the coarse domain, return
 704       if (i == 1) then
 705          get_nest_level = 1
 706          return
 707       end if
 708
 709       if (i > MAX_DOMAINS) then
 710          call mprintf(.true., ERROR, &
 711                       'get_nest_level() called with invalid grid ID of %i.',i1=i)
 712       end if
 713
 714       ! If not the coarse domain, then nesting level is at least 2
 715       ! Yes, this looks silly. But we do not have a grid_id array, so
 716       !    we must check on parent_id
 717       get_nest_level = 2
 718
 719       j = i
 720       do while (parent_id(j) /= 1)
 721          j = parent_id(j)
 722          get_nest_level = get_nest_level + 1
 723
 724          ! Sanity check
 725          if (get_nest_level > MAX_DOMAINS) then
 726             call mprintf(.true., ERROR, &
 727                          'Spooky nesting setup encountered in get_nest_level().')
 728          end if
 729       end do
 730
 731    end function get_nest_level
 732 #endif
 733
 734
 735    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 736    ! Name: select_domain
 737    !
 738    ! Purpose: This routine is used to select which nest x/y <-> lat/lon
 739    !   conversions will be with respect to. For example, selecting domain 2 will
 740    !   cause the llxy routine to compute x/y locations with respect to domain 2
 741    !   given a lat/lon.
 742    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 743    subroutine select_domain(domain_num)
 744
 745       implicit none
 746
 747       ! Arguments
 748       integer, intent(in) :: domain_num
 749
 750 #ifdef _GEOGRID
 751       if (domain_num > n_domains) then
 752          call mprintf(.true.,ERROR,'In select_domain(), selected domain is greater than n_domains.')
 753       end if
 754 #endif
 755 #ifdef _METGRID
 756       if (domain_num > 1) then
 757          call mprintf(.true.,ERROR,'In select_domain(), selected domain is greater than 1.')
 758       end if
 759 #endif
 760
 761       current_nest_number = domain_num
 762
 763    end subroutine select_domain
 764
 765
 766    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 767    ! Name: iget_selected_domain
 768    !
 769    ! Purpose: This function returns the number of the currently selected nest.
 770    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 771    function iget_selected_domain()
 772
 773       implicit none
 774
 775       ! Return value
 776       integer :: iget_selected_domain
 777
 778       iget_selected_domain = current_nest_number
 779
 780    end function iget_selected_domain
 781
 782
 783    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 784    ! Name: lltoxy
 785    !
 786    ! Purpose:
 787    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 788    subroutine lltoxy(xlat, xlon, x, y, stagger, comp_ll)
 789
 790       implicit none
 791
 792       ! Arguments
 793       integer, intent(in) :: stagger
 794       real, intent(in) :: xlat, xlon
 795       real, intent(out) :: x, y
 796       logical, optional, intent(in) :: comp_ll
 797
 798       ! Local variables
 799       logical :: save_comp_ll
 800
 801       ! Account for grid staggering
 802       if (stagger == HH) then
 803          proj_stack(current_nest_number)%stagger = HH
 804       else if (stagger == VV) then
 805          proj_stack(current_nest_number)%stagger = VV
 806       end if
 807
 808       if (present(comp_ll)) then
 809          save_comp_ll = proj_stack(current_nest_number)%comp_ll
 810          proj_stack(current_nest_number)%comp_ll = comp_ll
 811       end if
 812
 813       call latlon_to_ij(proj_stack(current_nest_number), xlat, xlon, x, y)
 814
 815       if (present(comp_ll)) then
 816          proj_stack(current_nest_number)%comp_ll = save_comp_ll
 817       end if
 818
 819       ! Account for grid staggering
 820       if (stagger == U) then
 821          x = x + 0.5
 822       else if (stagger == V) then
 823          y = y + 0.5
 824       else if (stagger == CORNER) then
 825          x = x + 0.5
 826          y = y + 0.5
 827       end if
 828
 829    end subroutine lltoxy
 830
 831
 832    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 833    ! Name: lltoxy
 834    !
 835    ! Purpose:
 836    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 837    subroutine xytoll(x, y, xlat, xlon, stagger, comp_ll)
 838
 839       implicit none
 840
 841       ! Arguments
 842       integer, intent(in) :: stagger
 843       real, intent(in) :: x, y
 844       real, intent(out) :: xlat, xlon
 845       logical, optional, intent(in) :: comp_ll
 846
 847       ! Local variables
 848       real :: rx, ry
 849       logical :: save_comp_ll
 850
 851       ! Account for grid staggering; we cannot modify x and y, so modify local
 852       !   copies of them
 853       if (stagger == U) then
 854          rx = x - 0.5
 855          ry = y
 856       else if (stagger == V) then
 857          rx = x
 858          ry = y - 0.5
 859       else if (stagger == HH) then
 860          proj_stack(current_nest_number)%stagger = HH
 861          rx = x
 862          ry = y
 863       else if (stagger == VV) then
 864          proj_stack(current_nest_number)%stagger = VV
 865          rx = x
 866          ry = y
 867       else if (stagger == CORNER) then
 868          proj_stack(current_nest_number)%stagger = CORNER
 869          rx = x - 0.5
 870          ry = y - 0.5
 871       else
 872          rx = x
 873          ry = y
 874       end if
 875
 876       if (present(comp_ll)) then
 877          save_comp_ll = proj_stack(current_nest_number)%comp_ll
 878          proj_stack(current_nest_number)%comp_ll = comp_ll
 879       end if
 880
 881       call ij_to_latlon(proj_stack(current_nest_number), rx, ry, xlat, xlon)
 882
 883       if (present(comp_ll)) then
 884          proj_stack(current_nest_number)%comp_ll = save_comp_ll
 885       end if
 886
 887    end subroutine xytoll
 888
 889 end module llxy_module