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