WPS/geogrid/util/plotter.F

   1 program plotter
   2
   3    implicit none
   4
   5    integer :: nx, ny
   6    integer :: i, j
   7    real :: start_r, start_g, start_b, end_r, end_g, end_b
   8    real :: lu, val, xlat, xlon, left, right, bottom, top, maxter, minter
   9
  10    call opngks
  11
  12    call gopwk(13, 41, 3)
  13
  14    call gscr(1, 0, 1.00, 1.00, 1.00)
  15    call gscr(1, 1, 0.00, 0.00, 0.00)
  16    call gscr(1, 2, 0.25, 0.25, 0.25)
  17    call gscr(1, 3, 1.00, 1.00, 0.50)
  18    call gscr(1, 4, 0.50, 1.00, 0.50)
  19    call gscr(1, 5, 1.00, 1.00, 0.00)
  20    call gscr(1, 6, 1.00, 1.00, 0.00)
  21    call gscr(1, 7, 0.50, 1.00, 0.50)
  22    call gscr(1, 8, 1.00, 1.00, 0.50)
  23    call gscr(1, 9, 0.50, 1.00, 0.50)
  24    call gscr(1,10, 0.50, 1.00, 0.50)
  25    call gscr(1,11, 1.00, 1.00, 0.50)
  26    call gscr(1,12, 0.00, 1.00, 0.00)
  27    call gscr(1,13, 0.00, 0.50, 0.00)
  28    call gscr(1,14, 0.00, 1.00, 0.00)
  29    call gscr(1,15, 0.00, 0.50, 0.00)
  30    call gscr(1,16, 0.00, 1.00, 0.00)
  31    call gscr(1,17, 0.50, 0.50, 1.00)
  32    call gscr(1,18, 0.00, 1.00, 0.00)
  33    call gscr(1,19, 0.00, 1.00, 0.00)
  34    call gscr(1,20, 0.75, 0.75, 0.75)
  35    call gscr(1,21, 0.75, 0.75, 0.75)
  36    call gscr(1,22, 0.00, 0.50, 0.00)
  37    call gscr(1,23, 0.75, 0.75, 0.75)
  38    call gscr(1,24, 0.75, 0.75, 0.75)
  39    call gscr(1,25, 1.00, 1.00, 1.00)
  40
  41    start_r = 0.00
  42    end_r   = 0.50
  43    start_g = 1.00
  44    end_g   = 0.25
  45    start_b = 0.00
  46    end_b   = 0.00
  47    do i=26,76
  48      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))
  49    end do
  50
  51    start_r = 0.50
  52    end_r   = 1.00
  53    start_g = 0.25
  54    end_g   = 1.00
  55    start_b = 0.00
  56    end_b   = 1.00
  57    do i=77,126
  58      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))
  59    end do
  60
  61    nx = 5324
  62    ny = 3344
  63
  64    left = 0.1
  65    right = 0.9
  66    bottom = 0.1
  67    top = 0.9
  68    call mappos(left,right,bottom,top)
  69    call mapstc('OU','US')
  70    call maproj('LC', 30., -98.00, 60.)
  71    call mapset('CO', 20.144764, -122.505325, 48.201309, -59.35916)
  72    call mapint()
  73
  74    open(42,file='lu.dat',form='formatted')
  75    open(43,file='lat.dat',form='formatted')
  76    open(44,file='lon.dat',form='formatted')
  77
  78    do j=1,ny
  79       do i=1,nx
  80          read(42,*) val
  81          read(43,*) xlat
  82          read(44,*) xlon
  83          call map_square(xlat, xlon, (right-left)/real(nx), (top-bottom)/real(nx), nint(val)+1)
  84       end do
  85    end do
  86
  87    close(42)
  88    close(43)
  89    close(44)
  90
  91    call maplot()
  92
  93    call frame()
  94
  95    open(41,file='lu.dat',form='formatted')
  96    open(42,file='ter.dat',form='formatted')
  97    open(43,file='lat.dat',form='formatted')
  98    open(44,file='lon.dat',form='formatted')
  99    maxter = -1000.
 100    minter = 10000.
 101    do j=1,ny
 102       do i=1,nx
 103          read(42,*) val
 104          if (val > maxter) maxter = val
 105          if (val < minter) minter = val
 106       end do
 107    end do
 108
 109    rewind(42)
 110
 111    do j=1,ny
 112       do i=1,nx
 113          read(41,*) lu
 114          read(42,*) val
 115          read(43,*) xlat
 116          read(44,*) xlon
 117          val = ((val-minter) * 99.)/(maxter-minter) + 26.
 118          if (nint(lu) ==  16) then
 119             call map_square(xlat, xlon, (right-left)/real(nx), (top-bottom)/real(nx), 17)
 120          else if (nint(lu) ==  1) then
 121             call map_square(xlat, xlon, (right-left)/real(nx), (top-bottom)/real(nx), 2)
 122          else
 123             call map_square(xlat, xlon, (right-left)/real(nx), (top-bottom)/real(nx), nint(val))
 124          end if
 125       end do
 126    end do
 127
 128    close(41)
 129    close(42)
 130    close(43)
 131    close(44)
 132
 133    call maplot()
 134
 135    call gclwk(13)
 136
 137    call clsgks
 138
 139 end program plotter
 140
 141
 142 subroutine map_square(rlat, rlon, width, height, colr)
 143
 144     implicit none
 145
 146     ! Arguments
 147     real :: rlat, rlon, width, height
 148     integer :: colr
 149
 150     ! Local variables
 151     real :: u, v
 152     real, dimension(4) :: xra, yra
 153     real, dimension(2000) :: dst
 154     integer, dimension(3000) :: ind
 155
 156     call maptrn(rlat, rlon, u, v)
 157
 158     u = u + (width/2.)
 159     v = v + (height/2.)
 160
 161     xra(1) = u-(width/2.)
 162     xra(2) = u+(width/2.)
 163     xra(3) = u+(width/2.)
 164     xra(4) = u-(width/2.)
 165
 166     yra(1) = v-(height/2.)
 167     yra(2) = v-(height/2.)
 168     yra(3) = v+(height/2.)
 169     yra(4) = v+(height/2.)
 170
 171     call sfsgfa(xra, yra, 4, dst, 2000, ind, 3000, colr)
 172
 173 end subroutine map_square