geogrid/src/gridinfo_module.F

   1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   2 ! MODULE GRIDINFO_MODULE
   3 !
   4 ! This module handles (i.e., acquires, stores, and makes available) all data
   5 !   describing the model domains to be processed.
   6 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   7 module gridinfo_module
   8
   9    use constants_module
  10    use misc_definitions_module
  11    use module_debug
  12
  13    ! Parameters
  14    integer, parameter :: MAX_DOMAINS = 21
  15
  16    ! Variables
  17    integer :: iproj_type, n_domains, io_form_output, dyn_opt
  18    integer, dimension(MAX_DOMAINS) :: parent_grid_ratio, parent_id, ixdim, jydim
  19    integer, dimension(MAX_DOMAINS) :: subgrid_ratio_x, subgrid_ratio_y
  20    real :: known_lat, known_lon, pole_lat, pole_lon, stand_lon, truelat1, truelat2, &
  21            known_x, known_y, dxkm, dykm, phi, lambda, ref_lat, ref_lon, ref_x, ref_y, &
  22            dlatdeg, dlondeg
  23    real, dimension(MAX_DOMAINS) :: parent_ll_x, parent_ll_y, parent_ur_x, parent_ur_y
  24    character (len=MAX_FILENAME_LEN) :: geog_data_path, opt_output_from_geogrid_path, opt_geogrid_tbl_path
  25
  26    character (len=128), dimension(MAX_DOMAINS) :: geog_data_res
  27    character (len=1) :: gridtype
  28    logical :: do_tiled_output
  29    logical, dimension(MAX_DOMAINS) :: grid_is_active
  30    integer :: debug_level
  31
  32    contains
  33
  34    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  35    ! Name: get_grid_params
  36    !
  37    ! Purpose: This subroutine retrieves all parameters regarding the model domains
  38    !    to be processed by geogrid.exe. This includes map parameters, domain
  39    !    size and location, and nest information.
  40    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  41    subroutine get_grid_params()
  42
  43       implicit none
  44
  45       ! Local variables
  46       integer :: i, j, max_dom, funit, io_form_geogrid, interval_seconds
  47       real :: dx, dy, rparent_gridpts
  48       integer, dimension(MAX_DOMAINS) :: i_parent_start, j_parent_start, &
  49                            s_we, e_we, s_sn, e_sn, &
  50                            start_year, start_month, start_day, start_hour, start_minute, start_second, &
  51                            end_year,   end_month,   end_day,   end_hour,   end_minute,   end_second
  52       character (len=128) :: map_proj
  53       character (len=128), dimension(MAX_DOMAINS) :: start_date, end_date
  54       character (len=3) :: wrf_core
  55       logical :: is_used, nest_outside
  56       logical, dimension(MAX_DOMAINS) :: active_grid
  57       logical :: nocolons
  58
  59       namelist /share/ wrf_core, max_dom, start_date, end_date, &
  60                         start_year, end_year, start_month, end_month, &
  61                         start_day, end_day, start_hour, end_hour, &
  62                         start_minute, end_minute, start_second, end_second, &
  63                         interval_seconds, &
  64                         io_form_geogrid, opt_output_from_geogrid_path, &
  65                         debug_level, active_grid, &
  66                         subgrid_ratio_x, subgrid_ratio_y, &
  67                         nocolons
  68       namelist /geogrid/ parent_id, parent_grid_ratio, &
  69                          i_parent_start, j_parent_start, s_we, e_we, s_sn, e_sn, &
  70                          map_proj, ref_x, ref_y, ref_lat, ref_lon, &
  71                          pole_lat, pole_lon, truelat1, truelat2, stand_lon, &
  72                          dx, dy, geog_data_res, geog_data_path, opt_geogrid_tbl_path
  73
  74       ! Set defaults for namelist variables
  75       debug_level = 0
  76       io_form_geogrid = 2
  77       wrf_core = 'ARW'
  78       max_dom = 1
  79       geog_data_path = 'NOT_SPECIFIED'
  80       ref_x = NAN
  81       ref_y = NAN
  82       ref_lat = NAN
  83       ref_lon = NAN
  84       dx = NAN
  85       dy = NAN
  86       map_proj = 'Lambert'
  87       pole_lat = 90.0
  88       pole_lon = 0.0
  89       truelat1 = NAN
  90       truelat2 = NAN
  91       stand_lon = NAN
  92       do i=1,MAX_DOMAINS
  93          geog_data_res(i) = 'default'
  94          parent_id(i) = 1
  95          parent_grid_ratio(i) = INVALID
  96          s_we(i) = 1
  97          e_we(i) = INVALID
  98          s_sn(i) = 1
  99          e_sn(i) = INVALID
 100          start_year(i) = 0
 101          start_month(i) = 0
 102          start_day(i) = 0
 103          start_hour(i) = 0
 104          start_minute(i) = 0
 105          start_second(i) = 0
 106          end_year(i) = 0
 107          end_month(i) = 0
 108          end_day(i) = 0
 109          end_hour(i) = 0
 110          end_minute(i) = 0
 111          end_second(i) = 0
 112          start_date(i) = '0000-00-00_00:00:00'
 113          end_date(i) = '0000-00-00_00:00:00'
 114          active_grid(i) = .true.
 115          subgrid_ratio_x(i) = 1
 116          subgrid_ratio_y(i) = 1
 117       end do
 118       opt_output_from_geogrid_path = './'
 119       opt_geogrid_tbl_path = 'geogrid/'
 120       interval_seconds = INVALID
 121       nocolons = .false.
 122
 123       ! Read parameters from Fortran namelist
 124       do funit=10,100
 125          inquire(unit=funit, opened=is_used)
 126          if (.not. is_used) exit
 127       end do
 128       open(funit,file='namelist.wps',status='old',form='formatted',err=1000)
 129       read(funit,share)
 130       read(funit,geogrid)
 131       close(funit)
 132
 133 ! BUG: More properly handle debug_level in module_debug
 134       if (debug_level.gt.100) then
 135          call set_debug_level(DEBUG)
 136       else
 137          call set_debug_level(WARN)
 138       end if
 139
 140       call mprintf(.true.,LOGFILE,'Using the following namelist variables:')
 141       call mprintf(.true.,LOGFILE,'&SHARE')
 142       call mprintf(.true.,LOGFILE,'  WRF_CORE         = %s',s1=wrf_core)
 143       call mprintf(.true.,LOGFILE,'  MAX_DOM          = %i',i1=max_dom)
 144       call mprintf(.true.,LOGFILE,'  START_YEAR       = %i',i1=start_year(1))
 145       do i=2,max_dom
 146          call mprintf(.true.,LOGFILE,'                   = %i',i1=start_year(i))
 147       end do
 148       call mprintf(.true.,LOGFILE,'  START_MONTH      = %i',i1=start_month(1))
 149       do i=2,max_dom
 150          call mprintf(.true.,LOGFILE,'                   = %i',i1=start_month(i))
 151       end do
 152       call mprintf(.true.,LOGFILE,'  START_DAY        = %i',i1=start_day(1))
 153       do i=2,max_dom
 154          call mprintf(.true.,LOGFILE,'                   = %i',i1=start_day(i))
 155       end do
 156       call mprintf(.true.,LOGFILE,'  START_HOUR       = %i',i1=start_hour(1))
 157       do i=2,max_dom
 158          call mprintf(.true.,LOGFILE,'                   = %i',i1=start_hour(i))
 159       end do
 160       call mprintf(.true.,LOGFILE,'  START_MINUTE     = %i',i1=start_minute(1))
 161       do i=2,max_dom
 162          call mprintf(.true.,LOGFILE,'                   = %i',i1=start_minute(i))
 163       end do
 164       call mprintf(.true.,LOGFILE,'  START_SECOND     = %i',i1=start_second(1))
 165       do i=2,max_dom
 166          call mprintf(.true.,LOGFILE,'                   = %i',i1=start_second(i))
 167       end do
 168       call mprintf(.true.,LOGFILE,'  END_YEAR         = %i',i1=end_year(1))
 169       do i=2,max_dom
 170          call mprintf(.true.,LOGFILE,'                   = %i',i1=end_year(i))
 171       end do
 172       call mprintf(.true.,LOGFILE,'  END_MONTH        = %i',i1=end_month(1))
 173       do i=2,max_dom
 174          call mprintf(.true.,LOGFILE,'                   = %i',i1=end_month(i))
 175       end do
 176       call mprintf(.true.,LOGFILE,'  END_DAY          = %i',i1=end_day(1))
 177       do i=2,max_dom
 178          call mprintf(.true.,LOGFILE,'                   = %i',i1=end_day(i))
 179       end do
 180       call mprintf(.true.,LOGFILE,'  END_HOUR         = %i',i1=end_hour(1))
 181       do i=2,max_dom
 182          call mprintf(.true.,LOGFILE,'                   = %i',i1=end_hour(i))
 183       end do
 184       call mprintf(.true.,LOGFILE,'  END_MINUTE       = %i',i1=end_minute(1))
 185       do i=2,max_dom
 186          call mprintf(.true.,LOGFILE,'                   = %i',i1=end_minute(i))
 187       end do
 188       call mprintf(.true.,LOGFILE,'  END_SECOND       = %i',i1=end_second(1))
 189       do i=2,max_dom
 190          call mprintf(.true.,LOGFILE,'                   = %i',i1=end_second(i))
 191       end do
 192       call mprintf(.true.,LOGFILE,'  START_DATE       = %s',s1=start_date(1))
 193       do i=2,max_dom
 194          call mprintf(.true.,LOGFILE,'                   = %s',s1=start_date(i))
 195       end do
 196       call mprintf(.true.,LOGFILE,'  END_DATE         = %s',s1=end_date(1))
 197       do i=2,max_dom
 198          call mprintf(.true.,LOGFILE,'                   = %s',s1=end_date(i))
 199       end do
 200       call mprintf(.true.,LOGFILE,'  INTERVAL_SECONDS = %i',i1=interval_seconds)
 201       call mprintf(.true.,LOGFILE,'  IO_FORM_GEOGRID  = %i',i1=io_form_geogrid)
 202       call mprintf(.true.,LOGFILE,'  OPT_OUTPUT_FROM_GEOGRID_PATH = %s',s1=opt_output_from_geogrid_path)
 203       call mprintf(.true.,LOGFILE,'  SUBGRID_RATIO_X  = %i',i1=subgrid_ratio_x(1))
 204       do i=2,max_dom
 205          call mprintf(.true.,LOGFILE,'                   = %i',i1=subgrid_ratio_x(i))
 206       enddo
 207       call mprintf(.true.,LOGFILE,'  SUBGRID_RATIO_Y  = %i',i1=subgrid_ratio_y(1))
 208       do i=2,max_dom
 209          call mprintf(.true.,LOGFILE,'                   = %i',i1=subgrid_ratio_y(i))
 210       enddo
 211
 212       call mprintf(.true.,LOGFILE,'  DEBUG_LEVEL      = %i',i1=debug_level)
 213       call mprintf(.true.,LOGFILE,'  ACTIVE_GRID      = %l',l1=active_grid(1))
 214       do i=2,max_dom
 215          call mprintf(.true.,LOGFILE,'                   = %l',l1=active_grid(i))
 216       end do
 217       call mprintf(.true.,LOGFILE,'  NOCOLONS         = %l',l1=nocolons)
 218       call mprintf(.true.,LOGFILE,'/')
 219
 220       call mprintf(.true.,LOGFILE,'&GEOGRID')
 221       call mprintf(.true.,LOGFILE,'  PARENT_ID         = %i',i1=parent_id(1))
 222       do i=2,max_dom
 223          call mprintf(.true.,LOGFILE,'                    = %i',i1=parent_id(i))
 224       end do
 225       call mprintf(.true.,LOGFILE,'  PARENT_GRID_RATIO = %i',i1=parent_grid_ratio(1))
 226       do i=2,max_dom
 227          call mprintf(.true.,LOGFILE,'                    = %i',i1=parent_grid_ratio(i))
 228       end do
 229       call mprintf(.true.,LOGFILE,'  I_PARENT_START    = %i',i1=i_parent_start(1))
 230       do i=2,max_dom
 231          call mprintf(.true.,LOGFILE,'                    = %i',i1=i_parent_start(i))
 232       end do
 233       call mprintf(.true.,LOGFILE,'  J_PARENT_START    = %i',i1=j_parent_start(1))
 234       do i=2,max_dom
 235          call mprintf(.true.,LOGFILE,'                    = %i',i1=j_parent_start(i))
 236       end do
 237       call mprintf(.true.,LOGFILE,'  S_WE              = %i',i1=s_we(1))
 238       do i=2,max_dom
 239          call mprintf(.true.,LOGFILE,'                    = %i',i1=s_we(i))
 240       end do
 241       call mprintf(.true.,LOGFILE,'  E_WE              = %i',i1=e_we(1))
 242       do i=2,max_dom
 243          call mprintf(.true.,LOGFILE,'                    = %i',i1=e_we(i))
 244       end do
 245       call mprintf(.true.,LOGFILE,'  S_SN              = %i',i1=s_sn(1))
 246       do i=2,max_dom
 247          call mprintf(.true.,LOGFILE,'                    = %i',i1=s_sn(i))
 248       end do
 249       call mprintf(.true.,LOGFILE,'  E_SN              = %i',i1=e_sn(1))
 250       do i=2,max_dom
 251          call mprintf(.true.,LOGFILE,'                    = %i',i1=e_sn(i))
 252       end do
 253       call mprintf(.true.,LOGFILE,'  GEOG_DATA_RES     = %s',s1=geog_data_res(1))
 254       do i=2,max_dom
 255          call mprintf(.true.,LOGFILE,'                    = %s',s1=geog_data_res(i))
 256       end do
 257       call mprintf(.true.,LOGFILE,'  DX                = %f',f1=dx)
 258       call mprintf(.true.,LOGFILE,'  DY                = %f',f1=dy)
 259       call mprintf(.true.,LOGFILE,'  MAP_PROJ          = %s',s1=map_proj)
 260       call mprintf(.true.,LOGFILE,'  POLE_LAT          = %f',f1=pole_lat)
 261       call mprintf(.true.,LOGFILE,'  POLE_LON          = %f',f1=pole_lon)
 262       call mprintf(.true.,LOGFILE,'  REF_LAT           = %f',f1=ref_lat)
 263       call mprintf(.true.,LOGFILE,'  REF_LON           = %f',f1=ref_lon)
 264       call mprintf(.true.,LOGFILE,'  REF_X             = %f',f1=ref_x)
 265       call mprintf(.true.,LOGFILE,'  REF_Y             = %f',f1=ref_y)
 266       call mprintf(.true.,LOGFILE,'  TRUELAT1          = %f',f1=truelat1)
 267       call mprintf(.true.,LOGFILE,'  TRUELAT2          = %f',f1=truelat2)
 268       call mprintf(.true.,LOGFILE,'  STAND_LON         = %f',f1=stand_lon)
 269       call mprintf(.true.,LOGFILE,'  GEOG_DATA_PATH    = %s',s1=geog_data_path)
 270       call mprintf(.true.,LOGFILE,'  OPT_GEOGRID_TBL_PATH = %s',s1=opt_geogrid_tbl_path)
 271       call mprintf(.true.,LOGFILE,'/')
 272
 273       dxkm = dx
 274       dykm = dy
 275
 276       known_lat = ref_lat
 277       known_lon = ref_lon
 278       known_x = ref_x
 279       known_y = ref_y
 280
 281       ! Convert wrf_core to uppercase letters
 282       do i=1,3
 283          if (ichar(wrf_core(i:i)) >= 97) wrf_core(i:i) = char(ichar(wrf_core(i:i))-32)
 284       end do
 285
 286       ! Before doing anything else, we must have a valid grid type
 287       gridtype = ' '
 288       if (wrf_core == 'ARW') then
 289          gridtype = 'C'
 290          dyn_opt = 2
 291       else if (wrf_core == 'NMM') then
 292          gridtype = 'E'
 293          dyn_opt = 4
 294       end if
 295
 296       ! Next, if this is NMM, we need to subtract 1 from the specified E_WE and E_SN;
 297       !    for some reason, these two variables need to be set to 1 larger than they
 298       !    really ought to be in the WRF namelist, so, to be consistent, we will do
 299       !    the same in the WPS namelist
 300       if (gridtype == 'E') then
 301          do i=1,max_dom
 302             e_we(i) = e_we(i) - 1
 303             e_sn(i) = e_sn(i) - 1
 304          end do
 305       end if
 306
 307       call mprintf(gridtype /= 'C' .and. gridtype /= 'E', ERROR, &
 308                    'A valid wrf_core must be specified in the namelist. '// &
 309                    'Currently, only "ARW" and "NMM" are supported.')
 310
 311       call mprintf(max_dom > MAX_DOMAINS, ERROR, &
 312                    'In namelist, max_dom must be <= %i. To run with more'// &
 313                    ' than %i domains, increase the MAX_DOMAINS parameter.', &
 314                    i1=MAX_DOMAINS, i2=MAX_DOMAINS)
 315
 316       ! Every domain must have a valid parent id
 317       do i=2,max_dom
 318          call mprintf(parent_id(i) <= 0 .or. parent_id(i) >= i, ERROR, &
 319                       'In namelist, the parent_id of domain %i must be in '// &
 320                       'the range 1 to %i.', i1=i, i2=i-1)
 321       end do
 322
 323       ! Check for valid geog_data_path
 324       j=1
 325       do i=1,len(geog_data_path)
 326          geog_data_path(j:j) = geog_data_path(i:i)
 327          if (geog_data_path(i:i) /= ' ') j = j + 1
 328       end do
 329       if (geog_data_path(1:1) == ' ') then
 330          call mprintf(.true.,ERROR,'In namelist, geog_data_path must be specified.')
 331       end if
 332       j = len_trim(geog_data_path)
 333       if (j >= MAX_FILENAME_LEN) then
 334          call mprintf(.true.,ERROR, &
 335                       'In namelist, geog_data_path must be strictly less '// &
 336                       'than 128 characters in length.')
 337       else
 338          if (geog_data_path(j:j) /= '/') then
 339             geog_data_path(j+1:j+1) = '/'
 340          end if
 341       end if
 342
 343       ! Paths need to end with a /
 344       j = len_trim(opt_geogrid_tbl_path)
 345       if (opt_geogrid_tbl_path(j:j) /= '/') then
 346          opt_geogrid_tbl_path(j+1:j+1) = '/'
 347       end if
 348
 349       j = len_trim(opt_output_from_geogrid_path)
 350       if (opt_output_from_geogrid_path(j:j) /= '/') then
 351          opt_output_from_geogrid_path(j+1:j+1) = '/'
 352       end if
 353
 354       ! Handle IOFORM+100 to do tiled IO
 355       if (io_form_geogrid > 100) then
 356          do_tiled_output = .true.
 357          io_form_geogrid = io_form_geogrid - 100
 358       else
 359          do_tiled_output = .false.
 360       end if
 361
 362       ! Check for valid io_form_geogrid
 363       if ( &
 364 #ifdef IO_BINARY
 365           io_form_geogrid /= BINARY .and. &
 366 #endif
 367 #ifdef IO_NETCDF
 368           io_form_geogrid /= NETCDF .and. &
 369 #endif
 370 #ifdef IO_GRIB1
 371           io_form_geogrid /= GRIB1 .and. &
 372 #endif
 373           .true. ) then
 374          call mprintf(.true.,WARN,'Valid io_form_geogrid values are:')
 375 #ifdef IO_BINARY
 376          call mprintf(.true.,WARN,'       %i (=BINARY)',i1=BINARY)
 377 #endif
 378 #ifdef IO_NETCDF
 379          call mprintf(.true.,WARN,'       %i (=NETCDF)',i1=NETCDF)
 380 #endif
 381 #ifdef IO_GRIB1
 382          call mprintf(.true.,WARN,'       %i (=GRIB1)',i1=GRIB1)
 383 #endif
 384          call mprintf(.true.,ERROR,'No valid value for io_form_geogrid was specified in the namelist.')
 385       end if
 386       io_form_output = io_form_geogrid
 387
 388       ! Convert map_proj to uppercase letters
 389       do i=1,len(map_proj)
 390          if (ichar(map_proj(i:i)) >= 97) map_proj(i:i) = char(ichar(map_proj(i:i))-32)
 391       end do
 392
 393       ! Assign parameters to module variables
 394       if ((index(map_proj, 'LAMBERT') /= 0) .and. &
 395           (len_trim(map_proj) == len('LAMBERT'))) then
 396          iproj_type = PROJ_LC
 397
 398       else if ((index(map_proj, 'MERCATOR') /= 0) .and. &
 399                (len_trim(map_proj) == len('MERCATOR'))) then
 400          iproj_type = PROJ_MERC
 401
 402       else if ((index(map_proj, 'POLAR') /= 0) .and. &
 403                (len_trim(map_proj) == len('POLAR'))) then
 404          iproj_type = PROJ_PS
 405
 406       else if ((index(map_proj, 'ROTATED_LL') /= 0) .and. &
 407                (len_trim(map_proj) == len('ROTATED_LL'))) then
 408          iproj_type = PROJ_ROTLL
 409
 410       else if ((index(map_proj, 'LAT-LON') /= 0) .and. &
 411                (len_trim(map_proj) == len('LAT-LON'))) then
 412          iproj_type = PROJ_CASSINI
 413
 414       else
 415          call mprintf(.true.,ERROR,&
 416                       'In namelist, invalid map_proj specified. Valid '// &
 417                       'projections are "lambert", "mercator", "polar", '// &
 418                       '"lat-lon", and "rotated_ll".')
 419       end if
 420
 421       ! For Cassini / lat-lon projections
 422       if (iproj_type == PROJ_CASSINI) then
 423
 424          ! If no dx,dy specified, assume global grid
 425          if (dx == NAN .and. dy == NAN) then
 426             dlondeg = 360. / (e_we(1)-s_we(1))   ! Here, we really do not want e_we-s_we+1
 427             dlatdeg = 180. / (e_sn(1)-s_sn(1))   ! Here, we really do not want e_we-s_we+1
 428             known_x = 1.
 429             known_y = 1.
 430             known_lon = stand_lon + dlondeg/2.
 431             known_lat = -90. + dlatdeg/2.
 432             dxkm = EARTH_RADIUS_M * PI * 2.0 / (e_we(1)-s_we(1))
 433             dykm = EARTH_RADIUS_M * PI       / (e_sn(1)-s_sn(1))
 434
 435          ! If dx,dy specified, however, assume regional grid
 436          else
 437             dlatdeg = dy
 438             dlondeg = dx
 439             dxkm = dlondeg * EARTH_RADIUS_M * PI * 2.0 / 360.0
 440             dykm = dlatdeg * EARTH_RADIUS_M * PI * 2.0 / 360.0
 441             if (known_lat == NAN .or. known_lon == NAN) then
 442                call mprintf(.true.,ERROR,'For lat-lon projection, if dx/dy are specified, '// &
 443                         'a regional domain is assumed, and a ref_lat,ref_lon must also be specified')
 444             end if
 445          end if
 446       end if
 447
 448       ! Manually set truelat2 = truelat1 if truelat2 not specified for Lambert
 449       if (iproj_type == PROJ_LC .and. truelat2 == NAN) then
 450          call mprintf ((truelat1 == NAN), ERROR, "No TRUELAT1 specified for Lambert conformal projection.")
 451          truelat2 = truelat1
 452       end if
 453
 454
 455       n_domains = max_dom
 456
 457       ! For C grid, let ixdim and jydim be the number of velocity points in
 458       !    each direction; for E grid, we will put the row and column back
 459       !    later; maybe this should be changed to be more clear, though.
 460       do i=1,n_domains
 461          ixdim(i) = e_we(i) - s_we(i) + 1
 462          jydim(i) = e_sn(i) - s_sn(i) + 1
 463       end do
 464
 465       if (gridtype == 'E') then
 466          phi = dykm*real(jydim(1)-1)/2.
 467          lambda = dxkm*real(ixdim(1)-1)
 468       end if
 469
 470       ! If the user hasn't supplied a known_x and known_y, assume the center of domain 1
 471       if (gridtype == 'E' .and. (known_x /= NAN .or. known_y /= NAN)) then
 472          call mprintf(.true.,WARN, &
 473                       'Namelist variables ref_x and ref_y cannot be used for NMM grids.'// &
 474                       ' (ref_lat, ref_lon) will refer to the center of the coarse grid.')
 475       else if (gridtype == 'C') then
 476          if (known_x == NAN .and. known_y == NAN) then
 477             known_x = ixdim(1) / 2.
 478             known_y = jydim(1) / 2.
 479          else if (known_x == NAN .or. known_y == NAN) then
 480             call mprintf(.true.,ERROR, &
 481                       'In namelist.wps, neither or both of ref_x, ref_y must be specified.')
 482          end if
 483       end if
 484
 485       ! Checks specific to E grid
 486       if (gridtype == 'E') then
 487
 488          ! E grid supports only the rotated lat/lon projection
 489          if (iproj_type /= PROJ_ROTLL) then
 490             call mprintf(.true., WARN, &
 491                          'For the NMM core, projection type must be rotated '// &
 492                          'lat/lon (map_proj=rotated_ll)')
 493             call mprintf(.true.,WARN,'Projection will be set to rotated_ll')
 494             iproj_type = PROJ_ROTLL
 495          end if
 496
 497          ! In the following check, add back the 1 that we had to subtract above
 498          !   for the sake of being consistent with WRF namelist
 499          call mprintf(mod(e_sn(1)+1,2) /= 0, ERROR, &
 500                       'For the NMM core, E_SN must be an even number for grid %i.', i1=1)
 501
 502          do i=2,n_domains
 503             call mprintf((parent_grid_ratio(i) /= 3), WARN, &
 504                          'For the NMM core, the parent_grid_ratio must be 3 for '// &
 505                          'domain %i. A ratio of 3 will be assumed.', i1=i)
 506             parent_grid_ratio(i) = 3
 507          end do
 508
 509          ! Check that nests have an acceptable number of grid points in each dimension
 510 !         do i=2,n_domains
 511 !            rparent_gridpts = real(ixdim(i)+2)/real(parent_grid_ratio(i))
 512 !            if (floor(rparent_gridpts) /= ceiling(rparent_gridpts)) then
 513 !               call mprintf(.true.,ERROR,'For nest %i, e_we must be 3n-2 '// &
 514 !                            'for some integer n > 1.', &
 515 !                            i1=i)
 516 !            end if
 517 !            rparent_gridpts = real(jydim(i)+2)/real(parent_grid_ratio(i))
 518 !            if (floor(rparent_gridpts) /= ceiling(rparent_gridpts)) then
 519 !               call mprintf(.true.,ERROR,'For nest %i, e_sn must be 3n-2 '// &
 520 !                            'for some odd integer n > 1.', &
 521 !                            i1=i)
 522 !            end if
 523 !         end do
 524
 525          do i=2,n_domains
 526             parent_ll_x(i) = 1.
 527             parent_ll_y(i) = 1.
 528          end do
 529
 530       ! Checks specific to C grid
 531       else if (gridtype == 'C') then
 532
 533          ! C grid does not support the rotated lat/lon projection
 534          call mprintf((iproj_type == PROJ_ROTLL), ERROR, &
 535                       'Rotated lat/lon projection is not supported for the ARW core. '// &
 536                       'Valid projecitons are "lambert", "mercator", "polar", and "lat-lon".')
 537
 538          ! Check that nests have an acceptable number of grid points in each dimension
 539          do i=2,n_domains
 540             rparent_gridpts = real(ixdim(i)-1)/real(parent_grid_ratio(i))
 541             if (floor(rparent_gridpts) /= ceiling(rparent_gridpts)) then
 542                call mprintf(.true.,ERROR,'For nest %i, (e_we-s_we+1) must be one greater '// &
 543                             'than an integer multiple of the parent_grid_ratio of %i.', &
 544                             i1=i, i2=parent_grid_ratio(i))
 545             end if
 546             rparent_gridpts = real(jydim(i)-1)/real(parent_grid_ratio(i))
 547             if (floor(rparent_gridpts) /= ceiling(rparent_gridpts)) then
 548                call mprintf(.true.,ERROR,'For nest %i, (e_sn-s_sn+1) must be one greater '// &
 549                             'than an integer multiple of the parent_grid_ratio of %i.', &
 550                             i1=i, i2=parent_grid_ratio(i))
 551             end if
 552          end do
 553
 554          ! Check that a nest does not extend outside of its parent grid
 555          nest_outside = .false.
 556          do i=2,n_domains
 557             if (i_parent_start(i) >= ixdim(parent_id(i))) then
 558                call mprintf(.true.,WARN,'Nest %i cannot have i_parent_start outside of parent domain.',i1=i)
 559                nest_outside = .true.
 560             else if (i_parent_start(i) + (ixdim(i) - 1)/parent_grid_ratio(i) > ixdim(parent_id(i))) then
 561                call mprintf(.true.,WARN,'Nest %i extends beyond its parent grid in the west-east direction.',i1=i)
 562                call mprintf(.true.,WARN,'   Maximum allowable e_we for current i/j_parent_start is %i.', &
 563                             i1=(ixdim(parent_id(i))-i_parent_start(i))*parent_grid_ratio(i)+1 )
 564                nest_outside = .true.
 565             end if
 566             if (j_parent_start(i) >= jydim(parent_id(i))) then
 567                call mprintf(.true.,WARN,'Nest %i cannot have j_parent_start outside of parent domain.',i1=i)
 568                nest_outside = .true.
 569             else if (j_parent_start(i) + (jydim(i) - 1)/parent_grid_ratio(i) > jydim(parent_id(i))) then
 570                call mprintf(.true.,WARN,'Nest %i extends beyond its parent grid in the south-north direction.',i1=i)
 571                call mprintf(.true.,WARN,'   Maximum allowable e_sn for current i/j_parent_start is %i.', &
 572                             i1=(jydim(parent_id(i))-j_parent_start(i))*parent_grid_ratio(i)+1 )
 573                nest_outside = .true.
 574             end if
 575          end do
 576          if (nest_outside) then
 577             call mprintf(.true.,ERROR,'One or more nested domains extend beyond their parent domains.')
 578          end if
 579
 580          do i=1,n_domains
 581             parent_ll_x(i) = real(i_parent_start(i))
 582             parent_ll_y(i) = real(j_parent_start(i))
 583             parent_ur_x(i) = real(i_parent_start(i))+real(ixdim(i))/real(parent_grid_ratio(i))-1.
 584             parent_ur_y(i) = real(j_parent_start(i))+real(jydim(i))/real(parent_grid_ratio(i))-1.
 585             grid_is_active(i) = active_grid(i)
 586          end do
 587
 588       end if
 589
 590       return
 591
 592  1000 call mprintf(.true.,ERROR,'Error opening file namelist.wps')
 593
 594    end subroutine get_grid_params
 595
 596 end module gridinfo_module