hwrf-baseline-20120103-1354/geogrid/util/plotgrid/src/plotgrid.F

   1 program plotgrid
   2
   3    use input_module
   4
   5    implicit none
   6
   7    external ulpr
   8
   9    integer :: n, i, j, nx, ny
  10    integer :: istatus, start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
  11               start_mem_k, end_mem_k, dyn_opt, &
  12               west_east_dim, south_north_dim, bottom_top_dim, map_proj, is_water, num_land_cat, &
  13               is_ice, is_urban, isoilwater, grid_id, parent_id, i_parent_start, j_parent_start, &
  14               i_parent_end, j_parent_end, parent_grid_ratio, &
  15               we_patch_s, we_patch_e, we_patch_s_stag, we_patch_e_stag, &
  16               sn_patch_s, sn_patch_e, sn_patch_s_stag, sn_patch_e_stag
  17    real :: width, height
  18    real :: dx, dy, cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, pole_lat, pole_lon
  19    real :: start_r, start_g, start_b, end_r, end_g, end_b
  20    real :: ll_lat, ll_lon, ur_lat, ur_lon
  21    real :: left, right, bottom, top, maxter, minter
  22    real :: rotang
  23    real, dimension(16) :: corner_lats, corner_lons
  24    real, dimension(10000) :: xcs, ycs
  25    integer, dimension(10) :: iai, iag
  26    integer, dimension(400000) :: iam
  27    integer, allocatable, dimension(:,:) :: lu
  28    real, allocatable, dimension(:,:) :: xlat, xlon, ter
  29    real, dimension(122000) :: rwrk
  30    real, pointer, dimension(:,:,:) :: real_array
  31    character (len=3) :: memorder
  32    character (len=25) :: crotang
  33    character (len=25) :: units
  34    character (len=46) :: desc
  35    character (len=128) :: init_date, cname, stagger, cunits, cdesc, title, startdate, grid_type, mminlu
  36    character (len=128), dimension(3) :: dimnames
  37
  38    call getarg(1,crotang)
  39    read (crotang,'(f)') rotang
  40
  41    write(6,*) 'Plotting with rotation angle ',rotang
  42
  43    call opngks
  44
  45    call gopwk(13, 41, 3)
  46
  47    call gscr(1, 0, 1.00, 1.00, 1.00)
  48    call gscr(1, 1, 0.00, 0.00, 0.00)
  49    call gscr(1, 2, 0.25, 0.25, 0.25)
  50    call gscr(1, 3, 1.00, 1.00, 0.50)
  51    call gscr(1, 4, 0.50, 1.00, 0.50)
  52    call gscr(1, 5, 1.00, 1.00, 0.00)
  53    call gscr(1, 6, 1.00, 1.00, 0.00)
  54    call gscr(1, 7, 0.50, 1.00, 0.50)
  55    call gscr(1, 8, 1.00, 1.00, 0.50)
  56    call gscr(1, 9, 0.50, 1.00, 0.50)
  57    call gscr(1,10, 0.50, 1.00, 0.50)
  58    call gscr(1,11, 1.00, 1.00, 0.50)
  59    call gscr(1,12, 0.00, 1.00, 0.00)
  60    call gscr(1,13, 0.00, 0.50, 0.00)
  61    call gscr(1,14, 0.00, 1.00, 0.00)
  62    call gscr(1,15, 0.00, 0.50, 0.00)
  63    call gscr(1,16, 0.00, 1.00, 0.00)
  64    call gscr(1,17, 0.50, 0.50, 1.00)
  65    call gscr(1,18, 0.00, 1.00, 0.00)
  66    call gscr(1,19, 0.00, 1.00, 0.00)
  67    call gscr(1,20, 0.75, 0.75, 0.75)
  68    call gscr(1,21, 0.75, 0.75, 0.75)
  69    call gscr(1,22, 0.00, 0.50, 0.00)
  70    call gscr(1,23, 0.75, 0.75, 0.75)
  71    call gscr(1,24, 0.75, 0.75, 0.75)
  72    call gscr(1,25, 1.00, 1.00, 1.00)
  73
  74    start_r = 0.00
  75    end_r   = 0.50
  76    start_g = 1.00
  77    end_g   = 0.25
  78    start_b = 0.00
  79    end_b   = 0.00
  80    do i=26,76
  81      call gscr(1,i,start_r+((end_r-start_r)/50.)*real(i-26),start_g+((end_g-start_g)/50.)*real(i-26),start_b+((end_b-start_b)/50.)*real(i-26))
  82    end do
  83
  84    start_r = 0.50
  85    end_r   = 1.00
  86    start_g = 0.25
  87    end_g   = 1.00
  88    start_b = 0.00
  89    end_b   = 1.00
  90    do i=77,126
  91      call gscr(1,i,start_r+((end_r-start_r)/50.)*real(i-77),start_g+((end_g-start_g)/50.)*real(i-77),start_b+((end_b-start_b)/50.)*real(i-77))
  92    end do
  93
  94    start_r = 0.80
  95    end_r   = 1.00
  96    start_g = 0.80
  97    end_g   = 1.00
  98    start_b = 0.80
  99    end_b   = 1.00
 100    do i=127,176
 101      call gscr(1,i,start_r+((end_r-start_r)/50.)*real(i-127),start_g+((end_g-start_g)/50.)*real(i-127),start_b+((end_b-start_b)/50.)*real(i-127))
 102    end do
 103
 104    call get_namelist_params()
 105
 106    do n=1,max_dom
 107       call input_init(n, istatus)
 108       if (istatus /= 0) then
 109          write(6,*) ' '
 110          write(6,*) 'Error: Could not open domain01 file.'
 111          write(6,*) ' '
 112          stop
 113       end if
 114
 115       call read_global_attrs(title, init_date, grid_type, dyn_opt, &
 116                              west_east_dim, south_north_dim, bottom_top_dim, &
 117                              we_patch_s, we_patch_e, we_patch_s_stag, we_patch_e_stag, &
 118                              sn_patch_s, sn_patch_e, sn_patch_s_stag, sn_patch_e_stag, &
 119                              map_proj, mminlu, num_land_cat, is_water, &
 120                              is_ice, is_urban, isoilwater, grid_id, parent_id, i_parent_start, j_parent_start, &
 121                              i_parent_end, j_parent_end, dx, dy, cen_lat, moad_cen_lat, cen_lon, &
 122                              stand_lon, truelat1, truelat2, pole_lat, pole_lon, parent_grid_ratio, &
 123                              corner_lats, corner_lons)
 124
 125       istatus = 0
 126       do while (istatus == 0)
 127         call read_next_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
 128                              start_mem_k, end_mem_k, cname, cunits, cdesc, &
 129                              memorder, stagger, dimnames, real_array, istatus)
 130         if (istatus == 0) then
 131
 132           if (index(cname, 'XLAT_M') /= 0) then
 133              nx = end_mem_i - start_mem_i + 1
 134              ny = end_mem_j - start_mem_j + 1
 135              allocate(xlat(nx,ny))
 136              xlat = real_array(:,:,1)
 137           else if (index(cname, 'XLONG_M') /= 0) then
 138              nx = end_mem_i - start_mem_i + 1
 139              ny = end_mem_j - start_mem_j + 1
 140              allocate(xlon(nx,ny))
 141              xlon = real_array(:,:,1)
 142           else if (index(cname, 'LU_INDEX') /= 0) then
 143              nx = end_mem_i - start_mem_i + 1
 144              ny = end_mem_j - start_mem_j + 1
 145              allocate(lu(nx,ny))
 146              lu = nint(real_array(:,:,1))
 147           end if
 148         end if
 149       end do
 150
 151       call input_close()
 152
 153       ll_lat = xlat(1,1)
 154       ll_lon = xlon(1,1)
 155       ur_lat = xlat(nx,ny)
 156       ur_lon = xlon(nx,ny)
 157 !      if (ur_lon < 0.) ur_lon = ur_lon + 360.0
 158
 159       if (n == 1) then
 160          left = 0.0
 161          right = 1.0
 162          bottom = 0.0
 163          top = 1.0
 164          call mappos(left,right,bottom,top)
 165
 166          call mapstc('OU','CO')
 167
 168          call maproj('CE', cen_lat, cen_lon, rotang)
 169 !         call maproj('LC', truelat1, stand_lon, truelat2)
 170 !         call maproj('ST', cen_lat, cen_lon, stand_lon)
 171          call mapset('CO', ll_lat, ll_lon, ur_lat, ur_lon)
 172          call mapint()
 173       end if
 174
 175       call mpsetr('GR', 10.0)
 176
 177       call maptrn(ll_lat, ll_lon, left, bottom)
 178       call maptrn(ur_lat, ur_lon, right, top)
 179
 180       width = 1.02*(right-left)/real(nx)
 181       height = 1.02*(top-bottom)/real(ny)
 182
 183       do j=1,ny
 184          do i=1,nx
 185             call map_square(xlat(i,j), xlon(i,j), width, height, lu(i,j)+1)
 186          end do
 187       end do
 188
 189       if (n > 1) then
 190          call gsplci(0)
 191          call lined(left-width/2.,bottom-height/2.,left-width/2.,top+height/2.)
 192          call lined(left-width/2.,top+height/2.,right+width/2.,top+height/2.)
 193          call lined(right+width/2.,top+height/2.,right+width/2.,bottom-height/2.)
 194          call lined(right+width/2.,bottom-height/2.,left-width/2.,bottom-height/2.)
 195          call sflush()
 196          call gsplci(1)
 197       end if
 198
 199       deallocate(xlat)
 200       deallocate(xlon)
 201       deallocate(lu)
 202    end do
 203
 204    call mplndr('Earth..3',4)
 205
 206    call arinam (iam,400000)
 207    call mapbla (iam)
 208    call arpram (iam,0,0,0)
 209
 210    call mapgrm (iam,xcs,ycs,10000,iai,iag,10,ulpr)
 211
 212    call frame()
 213
 214    do n=1,max_dom
 215       call input_init(n, istatus)
 216       if (istatus /= 0) then
 217          write(6,*) ' '
 218          write(6,*) 'Error: Could not open domain01 file.'
 219          write(6,*) ' '
 220          stop
 221       end if
 222
 223       call read_global_attrs(title, init_date, grid_type, dyn_opt, &
 224                              west_east_dim, south_north_dim, bottom_top_dim, &
 225                              we_patch_s, we_patch_e, we_patch_s_stag, we_patch_e_stag, &
 226                              sn_patch_s, sn_patch_e, sn_patch_s_stag, sn_patch_e_stag, &
 227                              map_proj, mminlu, num_land_cat, is_water, &
 228                              is_ice, is_urban, isoilwater, grid_id, parent_id, i_parent_start, j_parent_start, &
 229                              i_parent_end, j_parent_end, dx, dy, cen_lat, moad_cen_lat, cen_lon, &
 230                              stand_lon, truelat1, truelat2, pole_lat, pole_lon, parent_grid_ratio, &
 231                              corner_lats, corner_lons)
 232
 233       istatus = 0
 234       do while (istatus == 0)
 235         call read_next_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
 236                              start_mem_k, end_mem_k, cname, cunits, cdesc, &
 237                              memorder, stagger, dimnames, real_array, istatus)
 238         if (istatus == 0) then
 239
 240           if (index(cname, 'XLAT_M') /= 0) then
 241              nx = end_mem_i - start_mem_i + 1
 242              ny = end_mem_j - start_mem_j + 1
 243              allocate(xlat(nx,ny))
 244              xlat = real_array(:,:,1)
 245           else if (index(cname, 'XLONG_M') /= 0) then
 246              nx = end_mem_i - start_mem_i + 1
 247              ny = end_mem_j - start_mem_j + 1
 248              allocate(xlon(nx,ny))
 249              xlon = real_array(:,:,1)
 250           else if (index(cname, 'HGT_M') /= 0) then
 251              nx = end_mem_i - start_mem_i + 1
 252              ny = end_mem_j - start_mem_j + 1
 253              allocate(ter(nx,ny))
 254              ter = real_array(:,:,1)
 255           else if (index(cname, 'LU_INDEX') /= 0) then
 256              nx = end_mem_i - start_mem_i + 1
 257              ny = end_mem_j - start_mem_j + 1
 258              allocate(lu(nx,ny))
 259              lu = nint(real_array(:,:,1))
 260           end if
 261         end if
 262       end do
 263
 264       call input_close()
 265
 266       ll_lat = xlat(1,1)
 267       ll_lon = xlon(1,1)
 268       ur_lat = xlat(nx,ny)
 269       ur_lon = xlon(nx,ny)
 270
 271       if (n == 1) then
 272          left = 0.0
 273          right = 1.0
 274          bottom = 0.0
 275          top = 1.0
 276          call mappos(left,right,bottom,top)
 277
 278          call mapstc('OU','CO')
 279
 280          call maproj('CE', cen_lat, cen_lon, rotang)
 281 !         call maproj('LC', truelat1, stand_lon, truelat2)
 282 !         call maproj('ST', cen_lat, cen_lon, stand_lon)
 283          call mapset('CO', ll_lat, ll_lon, ur_lat, ur_lon)
 284          call mapint()
 285
 286          maxter = -10000.
 287          minter = 10000.
 288          do j=1,ny
 289             do i=1,nx
 290                if (ter(i,j) > maxter) maxter = ter(i,j)
 291                if (ter(i,j) < minter) minter = ter(i,j)
 292             end do
 293          end do
 294 !         maxter = 3348.42
 295       end if
 296
 297       call maptrn(ll_lat, ll_lon, left, bottom)
 298       call maptrn(ur_lat, ur_lon, right, top)
 299
 300       width = 1.02*(right-left)/real(nx)
 301       height = 1.02*(top-bottom)/real(ny)
 302
 303       do j=1,ny
 304          do i=1,nx
 305             if (lu(i,j) ==  16) then
 306                ter(i,j) = ((ter(i,j)-minter) * 99.)/(maxter-minter) + 26.
 307                call map_square(xlat(i,j), xlon(i,j), width, height, 17)
 308             else if (lu(i,j) ==  1) then
 309                ter(i,j) = ((ter(i,j)-minter) * 99.)/(maxter-minter) + 26.
 310                call map_square(xlat(i,j), xlon(i,j), width, height, 2)
 311             else if (lu(i,j) ==  24) then
 312                ter(i,j) = ((ter(i,j)-minter) * 50.)/(3500.0-minter) + 127.
 313                call map_square(xlat(i,j), xlon(i,j), width, height, nint(ter(i,j)))
 314             else
 315                ter(i,j) = ((ter(i,j)-minter) * 99.)/(maxter-minter) + 26.
 316                call map_square(xlat(i,j), xlon(i,j), width, height, nint(ter(i,j)))
 317             end if
 318          end do
 319       end do
 320
 321       if (n > 1) then
 322          call gsplci(0)
 323          call lined(left-width/2.,bottom-height/2.,left-width/2.,top+height/2.)
 324          call lined(left-width/2.,top+height/2.,right+width/2.,top+height/2.)
 325          call lined(right+width/2.,top+height/2.,right+width/2.,bottom-height/2.)
 326          call lined(right+width/2.,bottom-height/2.,left-width/2.,bottom-height/2.)
 327          call sflush()
 328          call gsplci(1)
 329       end if
 330
 331       deallocate(xlat)
 332       deallocate(xlon)
 333       deallocate(ter)
 334       deallocate(lu)
 335    end do
 336
 337    call mplndr('Earth..3',4)
 338
 339    call arinam (iam,400000)
 340    call mapbla (iam)
 341    call arpram (iam,0,0,0)
 342
 343    call mapgrm (iam,xcs,ycs,10000,iai,iag,10,ulpr)
 344
 345    call gclwk(13)
 346
 347    call clsgks
 348
 349
 350    stop
 351
 352 end program plotgrid
 353
 354
 355 subroutine map_square(rlat, rlon, width, height, colr)
 356
 357     implicit none
 358
 359     ! Arguments
 360     real :: rlat, rlon, width, height
 361     integer :: colr
 362
 363     ! Local variables
 364     real :: u, v
 365     real, dimension(4) :: xra, yra
 366     real, dimension(2000) :: dst
 367     integer, dimension(3000) :: ind
 368
 369     call maptrn(rlat, rlon, u, v)
 370
 371     xra(1) = u-(width/2.)
 372     xra(2) = u+(width/2.)
 373     xra(3) = u+(width/2.)
 374     xra(4) = u-(width/2.)
 375
 376     yra(1) = v-(height/2.)
 377     yra(2) = v-(height/2.)
 378     yra(3) = v+(height/2.)
 379     yra(4) = v+(height/2.)
 380
 381     call sfsgfa(xra, yra, 4, dst, 2000, ind, 3000, colr)
 382
 383 end subroutine map_square
 384
 385
 386 subroutine ulpr(xcs,ycs,ncs,iai,iag,nai)
 387
 388    implicit none
 389
 390    integer, external :: mapaci
 391
 392    integer :: ncs, nai
 393    integer, dimension(nai) :: iai, iag
 394    real, dimension(ncs) :: xcs, ycs
 395
 396    integer :: itm
 397
 398    if (iai(1) >= 0 .and.iai(2) >= 0) then
 399       itm = max0(iai(1),iai(2))
 400       if (mapaci(itm) == 1) then
 401          if (ncs.gt.150) print * , 'ulpr - ncs too big - ',ncs
 402          call gpl(ncs,xcs,ycs)
 403       end if
 404    end if
 405
 406 end subroutine ulpr