util/src/plotgrids.F

   1 program plotgrids
   2
   3    use map_utils
   4
   5    implicit none
   6
   7    ! Parameters
   8    integer, parameter :: MAX_DOMAINS = 21
   9
  10    ! Variables
  11    integer :: iproj_type, n_domains, io_form_output, dyn_opt
  12    integer :: i, j, max_dom, funit, io_form_geogrid
  13    integer :: interval_seconds
  14
  15    integer, dimension(MAX_DOMAINS) :: parent_grid_ratio, parent_id, ixdim, jydim
  16    integer, dimension(MAX_DOMAINS) :: i_parent_start, j_parent_start, &
  17                         s_we, e_we, s_sn, e_sn, &
  18                         start_year, start_month, start_day, start_hour, &
  19                         end_year,   end_month,   end_day,   end_hour
  20    logical, dimension(MAX_DOMAINS) :: active_grid
  21
  22    real :: known_lat, known_lon, stand_lon, truelat1, truelat2, known_x, known_y, &
  23            dxkm, dykm, ref_lat, ref_lon, ref_x, ref_y
  24    real :: dx, dy
  25    real :: ri, rj, rlats, rlons, rlate, rlone
  26    real :: polat , rot
  27    real :: rparent_gridpts
  28    real :: xa,xb,ya,yb,xxa,xxy,yya,yyb
  29    real :: xs, xe, ys, ye
  30    integer :: jproj, jgrid, jlts, iusout, idot, ier
  31    integer :: ltype , idom
  32
  33    real, dimension(MAX_DOMAINS) :: parent_ll_x, parent_ll_y, parent_ur_x, parent_ur_y
  34
  35    character (len=128) :: geog_data_path, opt_output_from_geogrid_path, opt_geogrid_tbl_path
  36    character (len=128), dimension(MAX_DOMAINS) :: geog_data_res
  37    character (len=128) :: map_proj
  38    character (len=128), dimension(MAX_DOMAINS) :: start_date, end_date
  39    character (len=3) :: wrf_core
  40    character (len=1) :: gridtype
  41
  42    logical :: do_tiled_output
  43    integer :: debug_level
  44    logical :: is_used
  45    logical :: nocolons
  46    integer :: isize
  47
  48    type (proj_info) :: map_projection
  49
  50    namelist /share/ wrf_core, max_dom, start_date, end_date, &
  51                      start_year, end_year, start_month, end_month, &
  52                      start_day, end_day, start_hour, end_hour, &
  53                      interval_seconds, io_form_geogrid, opt_output_from_geogrid_path, &
  54                      debug_level, active_grid, nocolons
  55    namelist /geogrid/ parent_id, parent_grid_ratio, &
  56                       i_parent_start, j_parent_start, s_we, e_we, s_sn, e_sn, &
  57                       map_proj, ref_x, ref_y, ref_lat, ref_lon, &
  58                       truelat1, truelat2, stand_lon, dx, dy, &
  59                       geog_data_res, geog_data_path, opt_geogrid_tbl_path
  60
  61    ! Set defaults for namelist variables
  62    debug_level = 0
  63    io_form_geogrid = 2
  64    wrf_core = 'ARW'
  65    max_dom = 1
  66    geog_data_path = 'NOT_SPECIFIED'
  67    ref_x = NAN
  68    ref_y = NAN
  69    ref_lat = NAN
  70    ref_lon = NAN
  71    dx = 10000.
  72    dy = 10000.
  73    map_proj = 'Lambert'
  74    truelat1 = NAN
  75    truelat2 = NAN
  76    stand_lon = NAN
  77    do i=1,MAX_DOMAINS
  78       geog_data_res(i) = 'default'
  79       parent_id(i) = 1
  80       parent_grid_ratio(i) = INVALID
  81       s_we(i) = 1
  82       e_we(i) = INVALID
  83       s_sn(i) = 1
  84       e_sn(i) = INVALID
  85       start_year(i) = 0
  86       start_month(i) = 0
  87       start_day(i) = 0
  88       start_hour(i) = 0
  89       end_year(i) = 0
  90       end_month(i) = 0
  91       end_day(i) = 0
  92       end_hour(i) = 0
  93       start_date(i) = '0000-00-00_00:00:00'
  94       end_date(i) = '0000-00-00_00:00:00'
  95    end do
  96    opt_output_from_geogrid_path = './'
  97    opt_geogrid_tbl_path = 'geogrid/'
  98    interval_seconds = INVALID
  99
 100    ! Read parameters from Fortran namelist
 101    do funit=10,100
 102       inquire(unit=funit, opened=is_used)
 103       if (.not. is_used) exit
 104    end do
 105    open(funit,file='namelist.wps',status='old',form='formatted',err=1000)
 106    read(funit,share)
 107    read(funit,geogrid)
 108    close(funit)
 109
 110    dxkm = dx
 111    dykm = dy
 112
 113    known_lat = ref_lat
 114    known_lon = ref_lon
 115    known_x = ref_x
 116    known_y = ref_y
 117
 118    ! Convert wrf_core to uppercase letters
 119    do i=1,3
 120       if (ichar(wrf_core(i:i)) >= 97) wrf_core(i:i) = char(ichar(wrf_core(i:i))-32)
 121    end do
 122
 123    ! Before doing anything else, we must have a valid grid type
 124    gridtype = ' '
 125    if (wrf_core == 'ARW') then
 126       gridtype = 'C'
 127       dyn_opt = 2
 128    else if (wrf_core == 'NMM') then
 129       gridtype = 'E'
 130       dyn_opt = 4
 131    end if
 132
 133    if (gridtype /= 'C' .and. gridtype /= 'E') then
 134       write(6,*) 'A valid wrf_core must be specified in the namelist. '// &
 135                  'Currently, only "ARW" and "NMM" are supported.'
 136       stop
 137    end if
 138
 139    if (max_dom > MAX_DOMAINS) then
 140       write(6,*) 'In namelist, max_dom must be <= ',MAX_DOMAINS,'. To run with more'// &
 141                 ' than ',MAX_DOMAINS,' domains, increase the MAX_DOMAINS parameter.'
 142       stop
 143    end if
 144
 145    ! Every domain must have a valid parent id
 146    do i=2,max_dom
 147       if (parent_id(i) <= 0 .or. parent_id(i) >= i) then
 148          write(6,*) 'In namelist, the parent_id of domain ',i,' must be in '// &
 149                    'the range 1 to ',i-1
 150           stop
 151       end if
 152    end do
 153
 154    ! Convert map_proj to uppercase letters
 155    do i=1,len(map_proj)
 156       if (ichar(map_proj(i:i)) >= 97) map_proj(i:i) = char(ichar(map_proj(i:i))-32)
 157    end do
 158
 159    ! Assign parameters to module variables
 160    if ((index(map_proj, 'LAMBERT') /= 0) .and. &
 161        (len_trim(map_proj) == len('LAMBERT'))) then
 162       iproj_type = PROJ_LC
 163       rot=truelat1
 164       polat=truelat2
 165       jproj = 3
 166
 167    else if ((index(map_proj, 'MERCATOR') /= 0) .and. &
 168             (len_trim(map_proj) == len('MERCATOR'))) then
 169       iproj_type = PROJ_MERC
 170       rot=0.
 171       polat=0.
 172       jproj = 9
 173
 174    else if ((index(map_proj, 'POLAR') /= 0) .and. &
 175             (len_trim(map_proj) == len('POLAR'))) then
 176       iproj_type = PROJ_PS
 177       rot=0.
 178       polat=SIGN(90., ref_lat)
 179       jproj = 1
 180
 181    else if ((index(map_proj, 'ROTATED_LL') /= 0) .and. &
 182             (len_trim(map_proj) == len('ROTATED_LL'))) then
 183       iproj_type = PROJ_ROTLL
 184
 185    else
 186          write(6,*) 'In namelist, invalid map_proj specified. Valid '// &
 187                     'projections are "lambert", "mercator", "polar", '// &
 188                     'and "rotated_ll".'
 189          stop
 190    end if
 191
 192    n_domains = max_dom
 193
 194    do i=1,n_domains
 195       ixdim(i) = e_we(i) - s_we(i) + 1
 196       jydim(i) = e_sn(i) - s_sn(i) + 1
 197    end do
 198
 199    ! If the user hasn't supplied a known_x and known_y, assume the center of domain 1
 200    if (known_x == NAN) known_x = ixdim(1) / 2.
 201    if (known_y == NAN) known_y = jydim(1) / 2.
 202
 203    ! Checks specific to C grid
 204    if (gridtype == 'C') then
 205
 206       ! C grid does not support the rotated lat/lon projection
 207       if (iproj_type == PROJ_ROTLL) then
 208          write(6,*) 'Rotated lat/lon projection is not supported for the ARW core. '// &
 209                     'Valid projecitons are "lambert", "mercator", and "polar".'
 210          stop
 211       end if
 212
 213       ! Check that nests have an acceptable number of grid points in each dimension
 214       do i=2,n_domains
 215          rparent_gridpts = real(ixdim(i)-1)/real(parent_grid_ratio(i))
 216          if (floor(rparent_gridpts) /= ceiling(rparent_gridpts)) then
 217             write(6,*) 'For nest ',i,' (e_we-s_we+1) must be one greater than an '// &
 218                        'integer multiple of the parent_grid_ratio.'
 219             write(6,*) 'Current values are s_we(i) = ',s_we(i),' e_we(i) = ',e_we(i)
 220             isize = nint(real(ixdim(i)-1)/real(parent_grid_ratio(i)))
 221             write(6,*) 'An e_we = ',isize * parent_grid_ratio(i) + 1,' might work'
 222             stop
 223          end if
 224          rparent_gridpts = real(jydim(i)-1)/real(parent_grid_ratio(i))
 225          if (floor(rparent_gridpts) /= ceiling(rparent_gridpts)) then
 226             write(6,*) 'For nest ',i,' (e_sn-s_sn+1) must be one greater than an '// &
 227                        'integer multiple of the parent_grid_ratio.'
 228             write(6,*) ' Current values are, s_sn(i) = ',s_sn(i),' e_sn(i) = ',e_sn(i)
 229             isize = nint(real(jydim(i)-1)/real(parent_grid_ratio(i)))
 230             write(6,*) 'An e_sn = ',isize * parent_grid_ratio(i) + 1,' might work'
 231             stop
 232          end if
 233       end do
 234    end if
 235
 236    ! Checks specific to E grid
 237    if (gridtype == 'E') then
 238
 239       ! E grid supports only the rotated lat/lon projection
 240       if (iproj_type /= PROJ_ROTLL) then
 241          write(6,*) 'Rotated lat/lon is the only supported projection for the NMM core.'
 242          stop
 243       end if
 244
 245       ! Check that the parent_grid_ratio is set to 3 for all nests
 246       do i=2,n_domains
 247          if (parent_grid_ratio(i) /= 3) then
 248             write(6,*) 'The parent_grid_ratio must be set to 3 for the NMM core.'
 249             stop
 250          end if
 251       end do
 252
 253       CALL plot_e_grid ( ref_lat , -1. * ref_lon , &
 254                          dy , dx, &
 255                          n_domains , &
 256                          e_we , e_sn , &
 257                          parent_id , parent_grid_ratio , &
 258                          i_parent_start , j_parent_start )
 259       stop
 260    end if
 261
 262    do i=1,n_domains
 263       parent_ll_x(i) = real(i_parent_start(i))
 264       parent_ll_y(i) = real(j_parent_start(i))
 265       parent_ur_x(i) = real(i_parent_start(i))+real(ixdim(i))/real(parent_grid_ratio(i))-1.
 266       parent_ur_y(i) = real(j_parent_start(i))+real(jydim(i))/real(parent_grid_ratio(i))-1.
 267    end do
 268
 269    call map_init(map_projection)
 270
 271    call map_set(iproj_type, map_projection, &
 272                 lat1=known_lat, &
 273                 lon1=known_lon, &
 274                 knowni=known_x, &
 275                 knownj=known_y, &
 276                 dx=dx, &
 277                 stdlon=stand_lon, &
 278                 truelat1=truelat1, &
 279                 truelat2=truelat2, &
 280                 ixdim=ixdim(1), &
 281                 jydim=jydim(1))
 282
 283    call ij_to_latlon(map_projection, 0.5, 0.5, rlats, rlons)
 284    call ij_to_latlon(map_projection, real(e_we(1))-0.5, real(e_sn(1))-0.5, rlate, rlone)
 285
 286    call opngks
 287
 288    ! Set some colors
 289    call gscr(1, 0, 1.00, 1.00, 1.00)
 290    call gscr(1, 1, 0.00, 0.00, 0.00)
 291
 292    ! Do not grind them with details
 293    jgrid=10
 294    jlts=-2
 295    iusout=1
 296    idot=0
 297
 298    call supmap(jproj,polat,stand_lon,rot,&
 299                rlats,rlons,rlate,rlone, &
 300                jlts,jgrid,iusout,idot,ier)
 301
 302    call setusv('LW',1000)
 303    call perim(e_we(1)-1,1,e_sn(1)-1,1)
 304    call getset(xa,xb,ya,yb,xxa,xxy,yya,yyb,ltype)
 305    call set   (xa,xb,ya,yb, &
 306          1.,real(e_we(1)),1.,real(e_sn(1)),ltype)
 307
 308    do idom = 2 , max_dom
 309       call getxy ( xs, xe, ys, ye, &
 310                    idom , max_dom , &
 311                    e_we , e_sn , &
 312                    parent_id , parent_grid_ratio , &
 313                    i_parent_start , j_parent_start )
 314       call line ( xs , ys , xe , ys )
 315       call line ( xe , ys , xe , ye )
 316       call line ( xe , ye , xs , ye )
 317       call line ( xs , ye , xs , ys )
 318    end do
 319
 320    call frame
 321
 322    write(6,*) ' '
 323    write(6,*) 'Creating plot in NCAR Graphics metafile...'
 324    write(6,*) ' '
 325
 326    call clsgks
 327
 328    write(6,*) ' *** Successful completion of program plotgrids.exe *** '
 329
 330
 331    stop
 332
 333 1000 write(6,*) 'Error opening namelist.wps'
 334    stop
 335
 336 end program plotgrids
 337
 338 subroutine getxy ( xs, xe, ys, ye, &
 339                    dom_id , num_domains , &
 340                    e_we , e_sn , &
 341                    parent_id , parent_grid_ratio , &
 342                    i_parent_start , j_parent_start )
 343
 344    implicit none
 345
 346    integer , intent(in) :: dom_id
 347    integer , intent(in) :: num_domains
 348    integer , intent(in) , dimension(num_domains):: e_we , e_sn , &
 349                                                    parent_id , parent_grid_ratio , &
 350                                                    i_parent_start , j_parent_start
 351    real , intent(out) :: xs, xe, ys, ye
 352
 353
 354    !  local vars
 355
 356    integer :: idom
 357
 358    xs = 0.
 359    xe = e_we(dom_id) -1
 360    ys = 0.
 361    ye = e_sn(dom_id) -1
 362
 363    idom = dom_id
 364    compute_xy : DO
 365
 366       xs = (i_parent_start(idom) + xs  -1 ) / &
 367            real(parent_grid_ratio(parent_id(idom)))
 368       xe = xe / real(parent_grid_ratio(idom))
 369
 370       ys = (j_parent_start(idom) + ys  -1 ) / &
 371            real(parent_grid_ratio(parent_id(idom)))
 372       ye = ye / real(parent_grid_ratio(idom))
 373
 374       idom = parent_id(idom)
 375       if ( idom .EQ. 1 ) then
 376          exit compute_xy
 377       end if
 378
 379    END DO compute_xy
 380
 381    xs = xs + 1
 382    xe = xs + xe
 383    ys = ys + 1
 384    ye = ys + ye
 385
 386 end subroutine getxy
 387
 388 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 389 !!!!!!  E GRID MAP INFO BELOW     !!!!!!!!!!!!!!
 390 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 391
 392 SUBROUTINE plot_e_grid ( rlat0d , rlon0d , dphd , dlmd, &
 393                          n_domains , &
 394                          e_we , e_sn , &
 395                          parent_id , parent_grid_ratio , &
 396                          i_parent_start , j_parent_start )
 397
 398 !  This routine generates a gmeta file of the area covered by an Arakawa e-grid.
 399 !  We assume that NCAR Graphics has not been called yet (and will be closed
 400 !  upon exit).  The required input fields are as from the WPS namelist file.
 401
 402    IMPLICIT NONE
 403
 404 !  15 April 2005 NCEP/EMC
 405 !    The Code and some instructions are provided by Tom BLACK to
 406 !    NCAR/DTC Meral Demirtas
 407
 408 !  4 May 2005  NCAR/DTC Meral DEMIRTAS
 409 !    - An include file (plot_inc) is added to get
 410 !    Domain size: IM,JM
 411 !    Central latitute and longnitute: RLAT0D,RLON0D
 412 !    Horizontal resolution: DPHD, DLMD
 413
 414 !  Feb 2007 NCAR/MMM
 415 !    Turn into f90
 416 !    Add implicit none
 417 !    Remove non-mapping portions
 418 !    Make part of WPS domain plotting utility
 419
 420 !  Dec 2008 NCAR/DTC
 421 !    Pass additional arguments to enable plotting of nests
 422
 423    !  Input map parameters for E grid.
 424
 425    REAL , INTENT(IN)    :: rlat0d , & ! latitude of grid center (degrees)
 426                            rlon0d     ! longitude of grid center (degrees, times -1)
 427
 428    REAL , INTENT(IN)    :: dphd , &   ! angular distance between rows (degrees)
 429                            dlmd       ! angular distance between adjacent H and V points (degrees)
 430
 431    INTEGER , INTENT(in) :: n_domains  ! number of domains
 432    INTEGER , INTENT(in) , DIMENSION(n_domains):: e_we , &
 433                                                  e_sn , &
 434                                                  parent_id , &
 435                                                  parent_grid_ratio , &
 436                                                  i_parent_start , &
 437                                                  j_parent_start
 438
 439    !  Some local vars
 440
 441    REAL :: rlat1d , &
 442            rlon1d
 443
 444    INTEGER :: im, &     ! number of H points in odd rows
 445               jm , &    ! number of rows
 446               ngpwe , &
 447               ngpsn , &
 448               ilowl , &
 449               jlowl
 450
 451    INTEGER :: imt , imtjm
 452    REAL :: latlft,lonlft,latrgt,lonrgt
 453
 454    im = e_we(1)-1
 455    jm = e_sn(1)-1
 456
 457    imt=2*im-1
 458    imtjm=imt*jm
 459    rlat1d=rlat0d
 460    rlon1d=rlon0d
 461    ngpwe=2*im-1
 462    ngpsn=jm
 463
 464    !  Get lat and lon of left and right points.
 465
 466    CALL corners ( rlat1d,rlon1d,im,jm,rlat0d,rlon0d,dphd,dlmd,&
 467                   ngpwe,ngpsn,ilowl,jlowl,latlft,lonlft,latrgt,lonrgt)
 468
 469    !  With corner points, make map background.
 470
 471    CALL mapbkg_egrid ( imt,jm,ilowl,jlowl,ngpwe,ngpsn,&
 472                        rlat0d,rlon0d,latlft,lonlft,latrgt,lonrgt,&
 473                        dlmd,dphd,&
 474                        n_domains,&
 475                        e_we,e_sn,&
 476                        parent_id,parent_grid_ratio,&
 477                        i_parent_start,j_parent_start)
 478
 479 END SUBROUTINE plot_e_grid
 480
 481 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 482
 483 SUBROUTINE corners ( glatd,glond,im,jm,tph0d,tlm0d,dphd,dlmd,&
 484                      ngpwe,ngpsn,ilowl,jlowl,glatl,glonl,glatr,glonr)
 485
 486    IMPLICIT NONE
 487
 488    REAL , INTENT(IN) :: glatd,glond,tph0d,tlm0d,dphd,dlmd,&
 489                         glatl,glonl,glatr,glonr
 490
 491    INTEGER , INTENT(IN) :: im,jm,ngpwe,ngpsn
 492    INTEGER , INTENT(OUT) :: ilowl,jlowl
 493
 494    !  Local vars
 495
 496    REAL , PARAMETER :: d2r = 1.74532925E-2 , r2D = 1./D2R
 497
 498    REAL :: glat , glon , dph , dlm , tph0 , tlm0
 499    REAL :: x , y , z , tlat , tlon , tlat1 , tlat2 , tlon1 , tlon2
 500    REAL :: row , col
 501    REAL :: dlm1 , dlm2 , d1 , d2 , d3 , d4 , dmin
 502
 503    INTEGER :: jmt , ii , jj , iuppr , juppr
 504    INTEGER :: nrow , ncol
 505
 506    jmt = jm/2+1
 507
 508    !  Convert from geodetic to transformed coordinates (degrees).
 509
 510    glat = glatd * d2r
 511    glon = glond * d2r
 512    dph = dphd * d2r
 513    dlm = dlmd * d2r
 514    tph0 = tph0d * d2r
 515    tlm0 = tlm0d * d2r
 516
 517    x = COS(tph0) * COS(glat) * COS(glon-tlm0)+SIN(tph0) * SIN(glat)
 518    y = -COS(glat) * SIN(glon-tlm0)
 519    z = COS(tph0) * SIN(glat)-SIN(tph0) * COS(glat) * COS(glon-tlm0)
 520    tlat = r2d * ATAN(z/SQRT(x*x + y*y))
 521    tlon = r2d * ATAN(y/x)
 522
 523    !  Find the real (non-integer) row and column of the input location on
 524    !  the filled e-grid.
 525
 526    row = tlat/dphd+jmt
 527    col = tlon/dlmd+im
 528    nrow = INT(row)
 529    ncol = INT(col)
 530    tlat = tlat * d2r
 531    tlon = tlon * d2r
 532
 533 !               E2     E3
 534 !
 535 !
 536 !                  X
 537 !               E1     E4
 538
 539    tlat1 = (nrow-jmt) * dph
 540    tlat2 = tlat1+dph
 541    tlon1 = (ncol-im) * dlm
 542    tlon2 = tlon1+dlm
 543
 544    dlm1 = tlon-tlon1
 545    dlm2 = tlon-tlon2
 546
 547    d1 = ACOS(COS(tlat) * COS(tlat1) * COS(dlm1)+SIN(tlat) * SIN(tlat1))
 548    d2 = ACOS(COS(tlat) * COS(tlat2) * COS(dlm1)+SIN(tlat) * SIN(tlat2))
 549    d3 = ACOS(COS(tlat) * COS(tlat2) * COS(dlm2)+SIN(tlat) * SIN(tlat2))
 550    d4 = ACOS(COS(tlat) * COS(tlat1) * COS(dlm2)+SIN(tlat) * SIN(tlat1))
 551
 552    dmin = MIN(d1,d2,d3,d4)
 553
 554    IF      ( ABS(dmin-d1) .LT. 1.e-6 ) THEN
 555      ii = ncol
 556      jj = nrow
 557    ELSE IF ( ABS(dmin-d2) .LT. 1.e-6 ) THEN
 558      ii = ncol
 559      jj = nrow+1
 560    ELSE IF ( ABS(dmin-d3) .LT. 1.e-6 ) THEN
 561      ii = ncol+1
 562      jj = nrow+1
 563    ELSE IF ( ABS(dmin-d4) .LT. 1.e-6 ) THEN
 564      ii = ncol+1
 565      jj = nrow
 566    END IF
 567
 568    !  Now find the i and j of the lower left corner of the desired grid
 569    !  region and of the upper right.
 570
 571    ilowl = ii-ngpwe/2
 572    jlowl = jj-ngpsn/2
 573    iuppr = ii+ngpwe/2
 574    juppr = jj+ngpsn/2
 575
 576    !  Find their geodetic coordinates.
 577
 578    CALL e2t2g(ilowl,jlowl,im,jm,tph0d,tlm0d,dphd,dlmd,glatl,glonl)
 579    CALL e2t2g(iuppr,juppr,im,jm,tph0d,tlm0d,dphd,dlmd,glatr,glonr)
 580
 581 END SUBROUTINE corners
 582
 583 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 584
 585 SUBROUTINE mapbkg_egrid ( imt,jm,ilowl,jlowl,ngpwe,ngpsn,&
 586                           rlat0d,rlon0d,glatl,glonl,glatr,glonr,&
 587                           dlmd,dphd,&
 588                           n_domains,&
 589                           e_we,e_sn,&
 590                           parent_id,parent_grid_ratio,&
 591                           i_parent_start , j_parent_start )
 592
 593 !  IMPLICIT NONE
 594
 595    INTEGER , INTENT(in) :: n_domains
 596    INTEGER , INTENT(in) , DIMENSION(n_domains):: e_we , e_sn , &
 597                                                  parent_id , parent_grid_ratio , &
 598                                                  i_parent_start , j_parent_start
 599
 600    !  Some local vars
 601
 602    CHARACTER (LEN=97) :: string
 603    INTEGER :: i
 604    REAL :: xs, xe, ys, ye
 605
 606    !  Yet more center lon messing around, hoo boy.
 607
 608 !  clonx=180.-rlon0d
 609    clonx=-rlon0d
 610
 611    !  Open up NCAR Graphics
 612
 613    CALL opngks
 614    CALL gopwk(8,9,3)
 615    CALL gsclip(0)
 616
 617    !  Make the background white, and the foreground black.
 618
 619    CALL gscr ( 1 , 0 , 1., 1., 1. )
 620    CALL gscr ( 1 , 1 , 0., 0., 0. )
 621
 622    !  Line width default thickness.
 623
 624    CALL setusv('LW',1000)
 625
 626    !  Make map outline a solid line, not dots.
 627
 628    CALL mapsti('MV',8)
 629    CALL mapsti('DO',0)
 630
 631    !  Map outlines are political and states.
 632
 633    CALL mapstc('OU','PS')
 634
 635    !  Cylindrical equidistant.
 636
 637    CALL maproj('CE',rlat0d,clonx,0.)
 638
 639    !  Specify corner points.
 640
 641    CALL mapset('CO',glatl,glonl,glatr,glonr)
 642
 643    !  Lat lon lines every 5 degrees.
 644
 645    CALL mapsti('GR',5)
 646
 647    !  Map takes up this much real estate.
 648
 649    CALL mappos( 0.05 , 0.95 , 0.05 , 0.95 )
 650
 651    !  Initialize and draw map.
 652
 653    CALL mapint
 654    CALL mapdrw
 655
 656    !  Line width twice default thickness.
 657
 658    CALL setusv('LW',2000)
 659
 660    !  Add approx grid point tick marks
 661
 662    CALL perim(((imt+3)/2)-1,1,jm-1,1)
 663
 664    !  Line width default thickness.
 665
 666    CALL setusv('LW',1000)
 667
 668    !  Put on a quicky description.
 669
 670    WRITE ( string , FMT = '("E-GRID E_WE = ",I4,", E_SN = ",I4 , &
 671                             &", DX = ",F6.4,", DY = ",F6.4 , &
 672                             &", REF_LAT = ",F8.3,", REF_LON = ",F8.3)') &
 673                             (imt+3)/2,jm+1,dlmd,dphd,rlat0d,-1.*rlon0d
 674    CALL getset(xa,xb,ya,yb,xxa,xxy,yya,yyb,ltype)
 675    CALL pchiqu(xxa,yya-(yyb-yya)/20.,string,8.,0.,-1.)
 676    CALL set   (xa,xb,ya,yb,&
 677                1.,real(e_we(1)),1.,real(e_sn(1)),ltype)
 678
 679    !  Line width twice default thickness.
 680
 681    CALL setusv('LW',2000)
 682
 683    !  Draw a box for each nest.
 684
 685    do i=2 , n_domains
 686       call getxy ( xs, xe, ys, ye, &
 687                    i , n_domains , &
 688                    e_we , e_sn , &
 689                    parent_id , parent_grid_ratio , &
 690                    i_parent_start , j_parent_start )
 691
 692       call line ( xs , ys , xe , ys )
 693       call line ( xe , ys , xe , ye )
 694       call line ( xe , ye , xs , ye )
 695       call line ( xs , ye , xs , ys )
 696    end do
 697
 698    CALL frame
 699
 700    !  Close workstation and NCAR Grpahics.
 701
 702    CALL gclwk(8)
 703    CALL clsgks
 704
 705 END SUBROUTINE mapbkg_egrid
 706
 707 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 708
 709 SUBROUTINE e2t2g ( ncol,nrow,im,jm,tph0d,tlm0d,dphd,dlmd,glatd,glond)
 710
 711 !  IMPLICIT NONE
 712
 713    REAL , PARAMETER :: D2R=1.74532925E-2 , R2D=1./D2R
 714
 715    DPH=DPHD*D2R
 716    DLM=DLMD*D2R
 717    TPH0=TPH0D*D2R
 718    TLM0=TLM0D*D2R
 719
 720 !***  FIND THE TRANSFORMED LAT (POSITIVE NORTH) AND LON (POSITIVE EAST)
 721
 722    TLATD=(NROW-(JM+1)/2)*DPHD
 723    TLOND=(NCOL-IM)*DLMD
 724
 725 !***  NOW CONVERT TO GEODETIC LAT (POSITIVE NORTH) AND LON (POSITIVE EAST)
 726
 727    TLATR=TLATD*D2R
 728    TLONR=TLOND*D2R
 729    ARG1=SIN(TLATR)*COS(TPH0)+COS(TLATR)*SIN(TPH0)*COS(TLONR)
 730    GLATR=ASIN(ARG1)
 731    GLATD=GLATR*R2D
 732    ARG2=COS(TLATR)*COS(TLONR)/(COS(GLATR)*COS(TPH0))- &
 733         TAN(GLATR)*TAN(TPH0)
 734    IF(ABS(ARG2).GT.1.)ARG2=ABS(ARG2)/ARG2
 735    FCTR=1.
 736    IF(TLOND.GT.0.)FCTR=-1.
 737    GLOND=TLM0D+FCTR*ACOS(ARG2)*R2D
 738    GLOND=-GLOND
 739
 740 END SUBROUTINE e2t2g