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