WPS/geogrid/util/plotgrid/src/plotgrid.F

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